source/cluster/hc.f90

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