!Hierarchical clustering using stored dissimilarities (Fortran) !hct.f, driver program !hc.f, routines !hcass.f, cluster assignments routine !hcden.f, draw part of dendrogram !To compile and link: f77 hct.f hc.f hcass.f hcden.f -o hct !To run on spectr.dat (hard-wired): hct !Hierarchical clustering without storage of dissimilarities (Fortran) !hcon2t.f, driver program !hcon2.f, routines !To compile and link: f77 hcon2t.f hcon2.f hcass.f hcden.f -o hcon2t !To run on spectr.dat (hard-wired): hcon2t ! REAL DATA(18,16),CRIT(18),MEMBR(18) ! REAL CRITVAL(9) ! INTEGER IA(18),IB(18) ! INTEGER ICLASS(18,9),HVALS(9) ! INTEGER IORDER(9),HEIGHT(9) ! DIMENSION NN(18),DISNN(18) ! REAL D(153) ! LOGICAL FLAG(18) Integer, parameter :: N=36, M=7, NN2=N*(N-1)/2, Lev=5 REAL DATA(N,M),CRIT(N),MEMBR(N) REAL CRITVAL(Lev) INTEGER IA(N),IB(N) INTEGER ICLASS(N,Lev),HVALS(Lev) INTEGER IORDER(Lev),HEIGHT(Lev) DIMENSION NN(N),DISNN(N) REAL D(NN2) LOGICAL FLAG(N) ! IN ABOVE, 18=N, 16=M, 9=LEV, 153=N(N-1)/2. ! ! OPEN(UNIT=21,STATUS='OLD',FILE='spectr.dat') ! ! ! N = 18 ! M = 16 DO I=1,N READ(21,100)(DATA(I,J),J=1,M) ENDDO 100 FORMAT(8F7.1) LEN = (N*(N-1))/2 IOPT=1 CALL HC(N,M,LEN,IOPT,DATA,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,D) !LEV = 9 CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT) !CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL) END !--------------------------------------------------------------- ! C ! HIERARCHICAL CLUSTERING using (user-specified) criterion. C ! C ! Parameters: C ! C ! DATA(N,M) input data matrix, C ! DISS(LEN) dissimilarities in lower half diagonal C ! storage; LEN = N.N-1/2, C ! IOPT clustering criterion to be used, C ! IA, IB, CRIT history of agglomerations; dimensions C ! N, first N-1 locations only used, C ! MEMBR, NN, DISNN vectors of length N, used to store C ! cluster cardinalities, current nearest C ! neighbour, and the dissimilarity assoc. C ! with the latter. C ! FLAG boolean indicator of agglomerable obj./ C ! clusters. C ! C ! F. Murtagh, ESA/ESO/STECF, Garching, February 1986. C ! C !------------------------------------------------------------C SUBROUTINE HC(N,M,LEN,IOPT,DATA,IA,IB,CRIT,MEMBR,NN,DISNN, FLAG,DISS) REAL DATA(N,M),MEMBR(N),DISS(LEN) INTEGER IA(N),IB(N) REAL CRIT(N) DIMENSION NN(N),DISNN(N) LOGICAL FLAG(N) REAL INF DATA INF/1.E+20/ ! ! Initializations ! DO I=1,N MEMBR(I)=1. FLAG(I)=.TRUE. ENDDO NCL=N ! ! Construct dissimilarity matrix ! DO I=1,N-1 DO J=I+1,N IND=IOFFSET(N,I,J) DISS(IND)=0. DO K=1,M DISS(IND)=DISS(IND)+(DATA(I,K)-DATA(J,K))**2 ENDDO IF (IOPT.EQ.1) DISS(IND)=DISS(IND)/2. ! (Above is done for the case of the min. var. method ! where merging criteria are defined in terms of variances ! rather than distances.) ENDDO ENDDO ! ! Carry out an agglomeration - first create list of NNs ! DO I=1,N-1 DMIN=INF DO J=I+1,N IND=IOFFSET(N,I,J) IF (DISS(IND).GE.DMIN) GOTO 500 DMIN=DISS(IND) JM=J 500 CONTINUE ENDDO NN(I)=JM DISNN(I)=DMIN ENDDO ! 400 CONTINUE ! Next, determine least diss. using list of NNs DMIN=INF DO I=1,N-1 IF (.NOT.FLAG(I)) GOTO 600 IF (DISNN(I).GE.DMIN) GOTO 600 DMIN=DISNN(I) IM=I JM=NN(I) 600 CONTINUE ENDDO NCL=NCL-1 ! ! This allows an agglomeration to be carried out. ! I2=MIN0(IM,JM) J2=MAX0(IM,JM) IA(N-NCL)=I2 IB(N-NCL)=J2 CRIT(N-NCL)=DMIN ! ! Update dissimilarities from new cluster. ! FLAG(J2)=.FALSE. DMIN=INF DO K=1,N IF (.NOT.FLAG(K)) GOTO 800 IF (K.EQ.I2) GOTO 800 X=MEMBR(I2)+MEMBR(J2)+MEMBR(K) IF (I2.LT.K) THEN IND1=IOFFSET(N,I2,K) ELSE IND1=IOFFSET(N,K,I2) ENDIF IF (J2.LT.K) THEN IND2=IOFFSET(N,J2,K) ELSE IND2=IOFFSET(N,K,J2) ENDIF IND3=IOFFSET(N,I2,J2) XX=DISS(IND3) ! ! WARD'S MINIMUM VARIANCE METHOD - IOPT=1. ! IF (IOPT.EQ.1) THEN DISS(IND1)=(MEMBR(I2)+MEMBR(K))*DISS(IND1)+ & (MEMBR(J2)+MEMBR(K))*DISS(IND2)- & MEMBR(K)*XX DISS(IND1)=DISS(IND1)/X ENDIF ! SINGLE LINK METHOD - IOPT=2. IF (IOPT.EQ.2) THEN DISS(IND1)=MIN(DISS(IND1),DISS(IND2)) ENDIF ! COMPLETE LINK METHOD - IOPT=3. IF (IOPT.EQ.3) THEN DISS(IND1)=MAX(DISS(IND1),DISS(IND2)) ENDIF ! AVERAGE LINK (OR GROUP AVERAGE) METHOD - IOPT=4. IF (IOPT.EQ.4) THEN DISS(IND1)=(MEMBR(I2)*DISS(IND1)+MEMBR(J2)*DISS(IND2))/ & (MEMBR(I2)+MEMBR(J2)) ENDIF ! MCQUITTY'S METHOD - IOPT=5. IF (IOPT.EQ.5) THEN DISS(IND1)=0.5*DISS(IND1)+0.5*DISS(IND2) ENDIF ! MEDIAN (GOWER'S) METHOD - IOPT=6. IF (IOPT.EQ.6) THEN DISS(IND1)=0.5*DISS(IND1)+0.5*DISS(IND2)-0.25*XX ENDIF ! CENTROID METHOD - IOPT=7. IF (IOPT.EQ.7) THEN DISS(IND1)=(MEMBR(I2)*DISS(IND1)+MEMBR(J2)*DISS(IND2)- & MEMBR(I2)*MEMBR(J2)*XX/(MEMBR(I2)+MEMBR(J2)))/ & (MEMBR(I2)+MEMBR(J2)) ENDIF IF (I2.GT.K) GOTO 800 IF (DISS(IND1).GE.DMIN) GOTO 800 DMIN=DISS(IND1) JJ=K 800 CONTINUE ENDDO MEMBR(I2)=MEMBR(I2)+MEMBR(J2) DISNN(I2)=DMIN NN(I2)=JJ ! Update list of NNs insofar as this is required. DO I=1,N-1 IF (.NOT.FLAG(I)) GOTO 900 IF (NN(I).EQ.I2) GOTO 850 IF (NN(I).EQ.J2) GOTO 850 GOTO 900 850 CONTINUE ! (Redetermine NN of I:) DMIN=INF DO J=I+1,N IND=IOFFSET(N,I,J) IF (.NOT.FLAG(J)) GOTO 870 IF (I.EQ.J) GOTO 870 IF (DISS(IND).GE.DMIN) GOTO 870 DMIN=DISS(IND) JJ=J 870 CONTINUE ENDDO NN(I)=JJ DISNN(I)=DMIN 900 CONTINUE ENDDO ! Repeat previous steps until N-1 agglomerations carried out. IF (NCL.GT.1) GOTO 400 RETURN END !--------------------------------------------- FUNCTION IOFFSET(N,I,J) ! Map row I and column J of upper half diagonal symmetric matrix ! onto vector. IOFFSET=J+(I-1)*N-(I*(I+1))/2 RETURN END !--------------------------------------------------------------- ! Given a HIERARCHIC CLUSTERING, described as a sequence of C ! agglomerations, derive the assignments into clusters for the C ! top LEV-1 levels of the hierarchy. C ! Prepare also the required data for representing the C ! dendrogram of this top part of the hierarchy. C ! C ! Parameters: C ! C ! IA, IB, CRIT: vectors of dimension N defining the agglomer- C ! ations. C ! LEV: number of clusters in largest partition. C ! HVALS: vector of dim. LEV, used internally only. C ! ICLASS: array of cluster assignments; dim. N by LEV. C ! IORDER, CRITVAL, HEIGHT: vectors describing the dendrogram, C ! all of dim. LEV. C ! F. Murtagh, ESA/ESO/STECF, Garching, February 1986. C ! HISTORY C ! Bounds bug fix, Oct. 1990, F. Murtagh. C ! Inserted line "IF (LOC.GT.LEV) GOTO 58" on line 48. This was C ! occassioned by incorrect termination of this loop when I C ! reached its (lower) extremity, i.e. N-LEV. Without the C ! /CHECK=(BOUNDS) option on VAX/VMS compilation, this inserted C ! statement was not necessary. C !---------------------------------------------------------------C SUBROUTINE HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT) INTEGER IA(N),IB(N),ICLASS(N,LEV),HVALS(LEV),IORDER(LEV),HEIGHT(LEV) REAL CRIT(N),CRITVAL(LEV) ! Pick out the clusters which the N objects belong to, ! at levels N-2, N-3, ... N-LEV+1 of the hierarchy. ! The clusters are identified by the lowest seq. no. of ! their members. ! There are 2, 3, ... LEV clusters, respectively, for the ! above levels of the hierarchy. HVALS(1)=1 HVALS(2)=IB(N-1) LOC=3 DO 59 I=N-2,N-LEV,-1 DO 52 J=1,LOC-1 IF (IA(I).EQ.HVALS(J)) GOTO 54 52 CONTINUE HVALS(LOC)=IA(I) LOC=LOC+1 54 CONTINUE DO 56 J=1,LOC-1 IF (IB(I).EQ.HVALS(J)) GOTO 58 56 CONTINUE IF (LOC.GT.LEV) GOTO 58 HVALS(LOC)=IB(I) LOC=LOC+1 58 CONTINUE 59 CONTINUE DO 400 LEVEL=N-LEV,N-2 DO 200 I=1,N ICL=I DO 100 ILEV=1,LEVEL 100 IF (IB(ILEV).EQ.ICL) ICL=IA(ILEV) NCL=N-LEVEL ICLASS(I,NCL-1)=ICL 200 CONTINUE 400 CONTINUE DO 120 I=1,N DO 120 J=1,LEV-1 DO 110 K=2,LEV IF (ICLASS(I,J).NE.HVALS(K)) GOTO 110 ICLASS(I,J)=K GOTO 120 110 CONTINUE 120 CONTINUE WRITE (6,450) 450 FORMAT(4X,' SEQ NOS 2CL 3CL 4CL 5CL 6CL 7CL 8CL 9CL') WRITE (6,470) 470 FORMAT(4X,' ------- --- --- --- --- --- --- --- --- ----') DO 500 I=1,N WRITE (6,600) I,(ICLASS(I,J),J=1,Lev) 600 FORMAT(I11,15I4) 500 CONTINUE ! Determine an ordering of the LEV clusters (at level LEV-1) ! for later representation of the dendrogram. ! These are stored in IORDER. ! Determine the associated ordering of the criterion values ! for the vertical lines in the dendrogram. ! The ordinal values of these criterion values may be used in ! preference, and these are stored in HEIGHT. ! Finally, note that the LEV clusters are renamed so that they ! have seq. nos. 1 to LEV. ! IORDER(1)=IA(N-1) IORDER(2)=IB(N-1) CRITVAL(1)=0.0 CRITVAL(2)=CRIT(N-1) HEIGHT(1)=LEV HEIGHT(2)=LEV-1 LOC=2 DO 700 I=N-2,N-LEV+1,-1 DO 650 J=1,LOC IF (IA(I).EQ.IORDER(J)) THEN ! Shift rightwards and insert IB(I) beside IORDER(J): DO 630 K=LOC+1,J+1,-1 IORDER(K)=IORDER(K-1) CRITVAL(K)=CRITVAL(K-1) HEIGHT(K)=HEIGHT(K-1) 630 CONTINUE IORDER(J+1)=IB(I) CRITVAL(J+1)=CRIT(I) HEIGHT(J+1)=I-(N-LEV) LOC=LOC+1 ENDIF 650 CONTINUE 700 CONTINUE DO 705 I=1,LEV DO 703 J=1,LEV IF (HVALS(I).EQ.IORDER(J)) THEN IORDER(J)=I GOTO 705 ENDIF 703 CONTINUE 705 CONTINUE RETURN END !--------------------------------------------------------------- ! Construct a DENDROGRAM of the top 8 levels of C ! a HIERARCHIC CLUSTERING. C ! C ! Parameters: C ! C ! IORDER, HEIGHT, CRITVAL: vectors of length LEV C ! defining the dendrogram. C ! These are: the ordering of objects C ! along the bottom of the dendrogram C ! (IORDER); the height of the vertical C ! above each object, in ordinal values C ! (HEIGHT); and in real values (CRITVAL).C ! C ! NOTE: these vectors MUST have been set up with C ! LEV = 9 in the prior call to routine C ! HCASS. ! C ! F. Murtagh, ESA/ESO/STECF, Garching, Feb. 1986.C ! C !-------------------------------------------------C SUBROUTINE HCDEN(LEV,IORDER,HEIGHT,CRITVAL) CHARACTER*80 LINE INTEGER IORDER(LEV),HEIGHT(LEV) REAL CRITVAL(LEV) INTEGER OUT(27,27) INTEGER UP,ACROSS,BLANK DATA UP,ACROSS,BLANK/'|','-',' '/ DO I=1,27 DO J=1,27 OUT(I,J)=BLANK ENDDO ENDDO DO I=3,27,3 I2=I/3 J2=28-3*HEIGHT(I2) DO J=27,J2,-1 OUT(J,I)=UP ENDDO DO K=I,3,-1 I3=INT((K+2)/3) IF ( (28-HEIGHT(I3)*3).LT.J2) GOTO 100 OUT(J2,K)=ACROSS ENDDO 100 CONTINUE ENDDO IC=3 DO I=1,27 IF (I.EQ.IC+1) THEN IDUM=IC/3 IDUM=9-IDUM DO L=1,9 IF (HEIGHT(L).EQ.IDUM) GOTO 190 ENDDO 190 IDUM=L WRITE(6,200) CRITVAL(IDUM),(OUT(I,J),J=1,27) IC=IC+3 ELSE LINE = ' ' WRITE(6,210) (OUT(I,J),J=1,27) ENDIF 200 FORMAT(1H ,8X,F12.2,4X,27A1) 210 FORMAT(1H ,24X,27A1) ENDDO WRITE(6,250) WRITE(6,220)(IORDER(J),J=1,9) WRITE(6,250) 220 FORMAT(1H ,24X,9I3) WRITE(6,230) 230 FORMAT(1H ,13X,'CRITERION CLUSTERS 1 TO 9') WRITE(6,240) 240 FORMAT(1H ,13X,'VALUES. (TOP 8 LEVELS OF HIERARCHY).') 250 FORMAT(/) RETURN END !--------------------------------------------------------------- !---------------------------------------------------------------