source/cluster/wham/src-M/hc.f

   1 C***********************  Contents  ****************************************
   2 C* Sample driver program, VAX-11 Fortran; **********************************
   3 C* HC: O(n^2) time, O(n^2) space hierarchical clustering, Fortran 77 *******
   4 C* HCASS: determine cluster-memberships, Fortran 77. ***********************
   5 C* HCDEN: draw upper part of dendrogram, VAX-11 Fortran. *******************
   6 C* Sample data set: last 36 lines. *****************************************
   7 C***************************************************************************
   8 C      REAL DATA(18,16),CRIT(18),MEMBR(18)
   9 C      REAL CRITVAL(9)
  10 C      INTEGER IA(18),IB(18)
  11 C      INTEGER ICLASS(18,9),HVALS(9)
  12 C      INTEGER IORDER(9),HEIGHT(9)
  13 C      DIMENSION NN(18),DISNN(18)
  14 C      REAL D(153)
  15 C      LOGICAL FLAG(18)
  16 C IN ABOVE, 18=N, 16=M, 9=LEV, 153=N(N-1)/2.
  17 C
  18 C
  19 C       OPEN(UNIT=21,STATUS='OLD',FILE='SPECTR.DAT')
  20 C
  21 C
  22 C      N = 18
  23 C      M = 16
  24 C      DO I=1,N
  25 C        READ(21,100)(DATA(I,J),J=1,M)
  26 C      ENDDO
  27 C 100  FORMAT(8F7.1)
  28 C
  29 C
  30 C      LEN = (N*(N-1))/2
  31 C      IOPT=1
  32 C      CALL HC(N,M,LEN,IOPT,DATA,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,D)
  33 C
  34 C
  35 C      LEV = 9
  36 C      CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
  37 C
  38 C
  39 C      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
  40 C
  41 C
  42 C      END
  43 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
  44 C                                                            C
  45 C  HIERARCHICAL CLUSTERING using (user-specified) criterion. C
  46 C                                                            C
  47 C  Parameters:                                               C
  48 C                                                            C
  49 Cremoved  DATA(N,M)         input data matrix,               C
  50 C  DISS(LEN)         dissimilarities in lower half diagonal  C
  51 C                    storage; LEN = N.N-1/2,                 C
  52 C  IOPT              clustering criterion to be used,        C
  53 C  IA, IB, CRIT      history of agglomerations; dimensions   C
  54 C                    N, first N-1 locations only used,       C
  55 C  MEMBR, NN, DISNN  vectors of length N, used to store      C
  56 C                    cluster cardinalities, current nearest  C
  57 C                    neighbour, and the dissimilarity assoc. C
  58 C                    with the latter.                        C
  59 C  FLAG              boolean indicator of agglomerable obj./ C
  60 C                    clusters.                               C
  61 C                                                            C
  62 C  F. Murtagh, ESA/ESO/STECF, Garching, February 1986.       C
  63 C                                                            C
  64 C------------------------------------------------------------C
  65       SUBROUTINE HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,
  66      X                FLAG,DISS)
  67       REAL MEMBR(N)
  68       REAL DISS(LEN)
  69       INTEGER IA(N),IB(N)
  70       REAL CRIT(N)
  71       DIMENSION NN(N),DISNN(N)
  72       LOGICAL FLAG(N)
  73       REAL INF
  74       DATA INF/1.E+20/
  75 C
  76 C  Initializations
  77 C
  78       DO I=1,N
  79          MEMBR(I)=1.
  80          FLAG(I)=.TRUE.
  81       ENDDO
  82       NCL=N
  83 C
  84 C  Construct dissimilarity matrix
  85 C
  86       DO I=1,N-1
  87          DO J=I+1,N
  88             IND=IOFFSET(N,I,J)
  89 cinput            DISS(IND)=0.
  90 cinput            DO K=1,M
  91 cinput               DISS(IND)=DISS(IND)+(DATA(I,K)-DATA(J,K))**2
  92 cinput            ENDDO
  93             IF (IOPT.EQ.1) DISS(IND)=DISS(IND)/2.
  94 C           (Above is done for the case of the min. var. method
  95 C            where merging criteria are defined in terms of variances
  96 C            rather than distances.)
  97           ENDDO
  98        ENDDO
  99 C
 100 C  Carry out an agglomeration - first create list of NNs
 101 C
 102       DO I=1,N-1
 103          DMIN=INF
 104          DO J=I+1,N
 105             IND=IOFFSET(N,I,J)
 106             IF (DISS(IND).GE.DMIN) GOTO 500
 107                DMIN=DISS(IND)
 108                JM=J
 109   500    CONTINUE
 110          ENDDO
 111          NN(I)=JM
 112          DISNN(I)=DMIN
 113       ENDDO
 114 C
 115   400 CONTINUE
 116 C     Next, determine least diss. using list of NNs
 117       DMIN=INF
 118       DO I=1,N-1
 119          IF (.NOT.FLAG(I)) GOTO 600
 120          IF (DISNN(I).GE.DMIN) GOTO 600
 121             DMIN=DISNN(I)
 122             IM=I
 123             JM=NN(I)
 124   600    CONTINUE
 125       ENDDO
 126       NCL=NCL-1
 127 C
 128 C  This allows an agglomeration to be carried out.
 129 C
 130       I2=MIN0(IM,JM)
 131       J2=MAX0(IM,JM)
 132       IA(N-NCL)=I2
 133       IB(N-NCL)=J2
 134       CRIT(N-NCL)=DMIN
 135 C
 136 C  Update dissimilarities from new cluster.
 137 C
 138       FLAG(J2)=.FALSE.
 139       DMIN=INF
 140       DO K=1,N
 141          IF (.NOT.FLAG(K)) GOTO 800
 142          IF (K.EQ.I2) GOTO 800
 143          X=MEMBR(I2)+MEMBR(J2)+MEMBR(K)
 144          IF (I2.LT.K) THEN
 145                            IND1=IOFFSET(N,I2,K)
 146                       ELSE
 147                            IND1=IOFFSET(N,K,I2)
 148          ENDIF
 149          IF (J2.LT.K) THEN
 150                            IND2=IOFFSET(N,J2,K)
 151                       ELSE
 152                            IND2=IOFFSET(N,K,J2)
 153          ENDIF
 154          IND3=IOFFSET(N,I2,J2)
 155          XX=DISS(IND3)
 156 C
 157 C  WARD'S MINIMUM VARIANCE METHOD - IOPT=1.
 158 C
 159          IF (IOPT.EQ.1) THEN
 160             DISS(IND1)=(MEMBR(I2)+MEMBR(K))*DISS(IND1)+
 161      X                 (MEMBR(J2)+MEMBR(K))*DISS(IND2)-
 162      X                 MEMBR(K)*XX
 163             DISS(IND1)=DISS(IND1)/X
 164          ENDIF
 165 C
 166 C  SINGLE LINK METHOD - IOPT=2.
 167 C
 168          IF (IOPT.EQ.2) THEN
 169             DISS(IND1)=MIN(DISS(IND1),DISS(IND2))
 170          ENDIF
 171 C
 172 C  COMPLETE LINK METHOD - IOPT=3.
 173 C
 174          IF (IOPT.EQ.3) THEN
 175             DISS(IND1)=MAX(DISS(IND1),DISS(IND2))
 176          ENDIF
 177 C
 178 C  AVERAGE LINK (OR GROUP AVERAGE) METHOD - IOPT=4.
 179 C
 180          IF (IOPT.EQ.4) THEN
 181             DISS(IND1)=(MEMBR(I2)*DISS(IND1)+MEMBR(J2)*DISS(IND2))/
 182      X                 (MEMBR(I2)+MEMBR(J2))
 183          ENDIF
 184 C
 185 C  MCQUITTY'S METHOD - IOPT=5.
 186 C
 187          IF (IOPT.EQ.5) THEN
 188             DISS(IND1)=0.5*DISS(IND1)+0.5*DISS(IND2)
 189          ENDIF
 190 C
 191 C  MEDIAN (GOWER'S) METHOD - IOPT=6.
 192 C
 193          IF (IOPT.EQ.6) THEN
 194             DISS(IND1)=0.5*DISS(IND1)+0.5*DISS(IND2)-0.25*XX
 195          ENDIF
 196 C
 197 C  CENTROID METHOD - IOPT=7.
 198 C
 199          IF (IOPT.EQ.7) THEN
 200             DISS(IND1)=(MEMBR(I2)*DISS(IND1)+MEMBR(J2)*DISS(IND2)-
 201      X          MEMBR(I2)*MEMBR(J2)*XX/(MEMBR(I2)+MEMBR(J2)))/
 202      X          (MEMBR(I2)+MEMBR(J2))
 203             ENDIF
 204 C
 205          IF (I2.GT.K) GOTO 800
 206          IF (DISS(IND1).GE.DMIN) GOTO 800
 207             DMIN=DISS(IND1)
 208             JJ=K
 209   800    CONTINUE
 210       ENDDO
 211       MEMBR(I2)=MEMBR(I2)+MEMBR(J2)
 212       DISNN(I2)=DMIN
 213       NN(I2)=JJ
 214 C
 215 C  Update list of NNs insofar as this is required.
 216 C
 217       DO I=1,N-1
 218          IF (.NOT.FLAG(I)) GOTO 900
 219          IF (NN(I).EQ.I2) GOTO 850
 220          IF (NN(I).EQ.J2) GOTO 850
 221          GOTO 900
 222   850    CONTINUE
 223 C        (Redetermine NN of I:)
 224          DMIN=INF
 225          DO J=I+1,N
 226             IND=IOFFSET(N,I,J)
 227             IF (.NOT.FLAG(J)) GOTO 870
 228             IF (I.EQ.J) GOTO 870
 229             IF (DISS(IND).GE.DMIN) GOTO 870
 230                DMIN=DISS(IND)
 231                JJ=J
 232   870       CONTINUE
 233          ENDDO
 234          NN(I)=JJ
 235          DISNN(I)=DMIN
 236   900    CONTINUE
 237       ENDDO
 238 C
 239 C  Repeat previous steps until N-1 agglomerations carried out.
 240 C
 241       IF (NCL.GT.1) GOTO 400
 242 C
 243 C
 244       RETURN
 245       END
 246 C
 247 C
 248       FUNCTION IOFFSET(N,I,J)
 249 C  Map row I and column J of upper half diagonal symmetric matrix
 250 C  onto vector.
 251       IOFFSET=J+(I-1)*N-(I*(I+1))/2
 252       RETURN
 253       END
 254 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
 255 C                                                               C
 256 C  Given a HIERARCHIC CLUSTERING, described as a sequence of    C
 257 C  agglomerations, derive the assignments into clusters for the C
 258 C  top LEV-1 levels of the hierarchy.                           C
 259 C  Prepare also the required data for representing the          C
 260 C  dendrogram of this top part of the hierarchy.                C
 261 C                                                               C
 262 C  Parameters:                                                  C
 263 C                                                               C
 264 C  IA, IB, CRIT: vectors of dimension N defining the agglomer-  C
 265 C                 ations.                                       C
 266 C  LEV:          number of clusters in largest partition.       C
 267 C  HVALS:        vector of dim. LEV, used internally only.      C
 268 C  ICLASS:       array of cluster assignments; dim. N by LEV.   C
 269 C  IORDER, CRITVAL, HEIGHT: vectors describing the dendrogram,  C
 270 C                all of dim. LEV.                               C
 271 C                                                               C
 272 C  F. Murtagh, ESA/ESO/STECF, Garching, February 1986.          C
 273 C                                                               C
 274 C HISTORY                                                       C
 275 C                                                               C
 276 C Bounds bug fix, Oct. 1990, F. Murtagh.                        C
 277 C Inserted line "IF (LOC.GT.LEV) GOTO 58" on line 48.  This was C
 278 C occassioned by incorrect termination of this loop when I      C
 279 C reached its (lower) extremity, i.e. N-LEV.  Without the       C
 280 C /CHECK=(BOUNDS) option on VAX/VMS compilation, this inserted  C
 281 C statement was not necessary.                                  C
 282 C---------------------------------------------------------------C
 283       SUBROUTINE HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,
 284      X        CRITVAL,HEIGHT)
 285       include 'sizesclu.dat'
 286       include 'COMMON.IOUNITS'
 287       integer ICLASS(maxconf,maxconf-1)
 288       INTEGER IA(N),IB(N),HVALS(LEV),IORDER(LEV),
 289      X        HEIGHT(LEV)
 290       REAL CRIT(N),CRITVAL(LEV)
 291 C
 292 C  Pick out the clusters which the N objects belong to,
 293 C  at levels N-2, N-3, ... N-LEV+1 of the hierarchy.
 294 C  The clusters are identified by the lowest seq. no. of
 295 C  their members.
 296 C  There are 2, 3, ... LEV clusters, respectively, for the
 297 C  above levels of the hierarchy.
 298 C
 299       HVALS(1)=1
 300       HVALS(2)=IB(N-1)
 301       LOC=3
 302       DO 59 I=N-2,N-LEV,-1
 303          DO 52 J=1,LOC-1
 304             IF (IA(I).EQ.HVALS(J)) GOTO 54
 305   52     CONTINUE
 306          HVALS(LOC)=IA(I)
 307          LOC=LOC+1
 308   54     CONTINUE
 309          DO 56 J=1,LOC-1
 310             IF (IB(I).EQ.HVALS(J)) GOTO 58
 311   56     CONTINUE
 312          IF (LOC.GT.LEV) GOTO 58
 313          HVALS(LOC)=IB(I)
 314          LOC=LOC+1
 315   58     CONTINUE
 316   59  CONTINUE
 317 C
 318       DO 400 LEVEL=N-LEV,N-2
 319          DO 200 I=1,N
 320             ICL=I
 321             DO 100 ILEV=1,LEVEL
 322   100       IF (IB(ILEV).EQ.ICL) ICL=IA(ILEV)
 323             NCL=N-LEVEL
 324             ICLASS(I,NCL-1)=ICL
 325   200    CONTINUE
 326   400  CONTINUE
 327 C
 328       DO 120 I=1,N
 329       DO 120 J=1,LEV-1
 330       DO 110 K=2,LEV
 331       IF (ICLASS(I,J).NE.HVALS(K)) GOTO 110
 332          ICLASS(I,J)=K
 333          GOTO 120
 334   110 CONTINUE
 335   120 CONTINUE
 336 C
 337       WRITE (iout,450) (j,j=2,LEV)
 338   450 FORMAT(4X,' SEQ NOS',8(i2,'CL'),10000(i3,'CL'))
 339       WRITE (iout,470) (' ---',j=2,LEV)
 340   470 FORMAT(4X,' -------',10000a4)
 341       DO 500 I=1,N
 342       WRITE (iout,600) I,(ICLASS(I,J),J=1,LEV-1)
 343   600 FORMAT(I11,8I4,10000i5)
 344   500 CONTINUE
 345 C
 346 C  Determine an ordering of the LEV clusters (at level LEV-1)
 347 C  for later representation of the dendrogram.
 348 C  These are stored in IORDER.
 349 C  Determine the associated ordering of the criterion values
 350 C  for the vertical lines in the dendrogram.
 351 C  The ordinal values of these criterion values may be used in
 352 C  preference, and these are stored in HEIGHT.
 353 C  Finally, note that the LEV clusters are renamed so that they
 354 C  have seq. nos. 1 to LEV.
 355 C
 356       IORDER(1)=IA(N-1)
 357       IORDER(2)=IB(N-1)
 358       CRITVAL(1)=0.0
 359       CRITVAL(2)=CRIT(N-1)
 360       HEIGHT(1)=LEV
 361       HEIGHT(2)=LEV-1
 362       LOC=2
 363       DO 700 I=N-2,N-LEV+1,-1
 364          DO 650 J=1,LOC
 365             IF (IA(I).EQ.IORDER(J)) THEN
 366 C              Shift rightwards and insert IB(I) beside IORDER(J):
 367                DO 630 K=LOC+1,J+1,-1
 368                   IORDER(K)=IORDER(K-1)
 369                   CRITVAL(K)=CRITVAL(K-1)
 370                   HEIGHT(K)=HEIGHT(K-1)
 371   630          CONTINUE
 372                IORDER(J+1)=IB(I)
 373                 CRITVAL(J+1)=CRIT(I)
 374                 HEIGHT(J+1)=I-(N-LEV)
 375                LOC=LOC+1
 376             ENDIF
 377   650   CONTINUE
 378   700 CONTINUE
 379       DO 705 I=1,LEV
 380          DO 703 J=1,LEV
 381             IF (HVALS(I).EQ.IORDER(J)) THEN
 382                IORDER(J)=I
 383                GOTO 705
 384             ENDIF
 385   703    CONTINUE
 386   705 CONTINUE
 387 C
 388       RETURN
 389       END
 390 C+++++++++++++++++++++++++++++++++++++++++++++++++C
 391 C                                                 C
 392 C  Construct a DENDROGRAM of the top 8 levels of  C
 393 C  a HIERARCHIC CLUSTERING.                       C
 394 C                                                 C
 395 C  Parameters:                                    C
 396 C                                                 C
 397 C  IORDER, HEIGHT, CRITVAL: vectors of length LEV C
 398 C          defining the dendrogram.               C
 399 C          These are: the ordering of objects     C
 400 C          along the bottom of the dendrogram     C
 401 C          (IORDER); the height of the vertical   C
 402 C          above each object, in ordinal values   C
 403 C          (HEIGHT); and in real values (CRITVAL).C
 404 C                                                 C
 405 C  NOTE: these vectors MUST have been set up with C
 406 C        LEV = 9 in the prior call to routine     C
 407 C        HCASS.
 408 C                                                 C
 409 C  F. Murtagh, ESA/ESO/STECF, Garching, Feb. 1986.C
 410 C                                                 C
 411 C-------------------------------------------------C
 412       SUBROUTINE HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 413       include 'COMMON.IOUNITS'
 414       CHARACTER*80 LINE
 415       INTEGER IORDER(LEV),HEIGHT(LEV)
 416       REAL CRITVAL(LEV)
 417 c      INTEGER OUT(3*LEV,3*LEV)
 418 c      INTEGER UP,ACROSS,BLANK
 419       CHARACTER*1 OUT(3*LEV,3*LEV)
 420       CHARACTER*1 UP,ACROSS,BLANK
 421       DATA UP,ACROSS,BLANK/'|','-',' '/
 422 C
 423 C
 424       DO I=1,3*LEV
 425         DO J=1,3*LEV
 426           OUT(I,J)=BLANK
 427         ENDDO
 428       ENDDO
 429 C
 430 C
 431       DO I=3,3*LEV,3
 432          I2=I/3
 433 C
 434          J2=3*LEV+1-3*HEIGHT(I2)
 435          DO J=3*LEV,J2,-1
 436             OUT(J,I)=UP
 437          ENDDO
 438 C
 439          DO K=I,3,-1
 440             I3=INT((K+2)/3)
 441             IF ( (3*LEV+1-HEIGHT(I3)*3).LT.J2) GOTO 100
 442             OUT(J2,K)=ACROSS
 443          ENDDO
 444   100    CONTINUE
 445 C
 446       ENDDO
 447 C
 448 C
 449       IC=3
 450       DO I=1,3*LEV
 451       IF (I.EQ.IC+1) THEN
 452                    IDUM=IC/3
 453                    IDUM=LEV-IDUM
 454                    DO L=1,LEV
 455                       IF (HEIGHT(L).EQ.IDUM) GOTO 190
 456                    ENDDO
 457   190              IDUM=L
 458                    WRITE(iout,200) CRITVAL(IDUM),(OUT(I,J),J=1,3*LEV)
 459                    IC=IC+3
 460                    ELSE
 461                    LINE = ' '
 462                    WRITE(iout,210) (OUT(I,J),J=1,3*LEV)
 463       ENDIF
 464   200 FORMAT(1H ,8X,F12.2,4X,27000A1)
 465   210 FORMAT(1H ,24X,27000A1)
 466       ENDDO
 467       WRITE(iout,250)
 468       WRITE(iout,220)(IORDER(J),J=1,LEV)
 469       WRITE(iout,250)
 470   220 FORMAT(1H ,24X,9000I3)
 471       WRITE(iout,230) LEV
 472   230 FORMAT(1H ,13X,'CRITERION        CLUSTERS 1 TO ',i3)
 473       WRITE(iout,240) LEV-1
 474   240 FORMAT(1H ,13X,'VALUES.      (TOP ',i3,' LEVELS OF HIERARCHY).')
 475   250 FORMAT(/)
 476 C
 477 C
 478       RETURN
 479       END