--- /dev/null
+!*********************** Contents ****************************************
+!* Sample driver program, VAX-11 Fortran; **********************************
+!* HC: O(n^2) time, O(n^2) space hierarchical clustering, Fortran 77 *******
+!* HCASS: determine cluster-memberships, Fortran 77. ***********************
+!* HCDEN: draw upper part of dendrogram, VAX-11 Fortran. *******************
+!* Sample data set: last 36 lines. *****************************************
+!***************************************************************************
+! 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)
+! 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
+! C
+! HIERARCHICAL CLUSTERING using (user-specified) criterion. C
+! C
+! Parameters: C
+! C
+!removed 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
+ module hc_
+!-----------------------------------------------------------------------------
+ use io_units
+ use names
+ use clust_data
+ implicit none
+!-----------------------------------------------------------------------------
+!
+!
+!-----------------------------------------------------------------------------
+ contains
+!-----------------------------------------------------------------------------
+! hc.f
+!-----------------------------------------------------------------------------
+
+ SUBROUTINE HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,&
+ FLAG,DISS)
+ integer :: N,M,LEN,IOPT
+ REAL(kind=4) :: MEMBR(N)
+ REAL(kind=4) :: DISS(LEN)
+ INTEGER :: IA(N),IB(N)
+ REAL(kind=4) :: CRIT(N)
+ integer,DIMENSION(N) :: NN
+ real(kind=4),dimension(N) ::DISNN
+ LOGICAL :: FLAG(N)
+ REAL(kind=4) INF
+ DATA INF /1.E+20/
+ integer :: I,J,NCL,IND,IM,JM,I2,J2,K,IND1,IND2,IND3,JJ
+ real(kind=8) :: DMIN,X,XX
+!
+! 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)
+!input DISS(IND)=0.
+!input DO K=1,M
+!input DISS(IND)=DISS(IND)+(DATA(I,K)-DATA(J,K))**2
+!input 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 SUBROUTINE HC
+!-----------------------------------------------------------------------------
+!
+!
+ integer FUNCTION IOFFSET(N,I,J)
+! Map row I and column J of upper half diagonal symmetric matrix
+! onto vector.
+ integer :: N,I,J
+ IOFFSET=J+(I-1)*N-(I*(I+1))/2
+ RETURN
+ END FUNCTION IOFFSET
+!-----------------------------------------------------------------------------
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
+! C
+! 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
+! C
+! F. Murtagh, ESA/ESO/STECF, Garching, February 1986. C
+! C
+! HISTORY C
+! 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)
+! include 'sizesclu.dat'
+! include 'COMMON.IOUNITS'
+ integer :: N,LEV
+!el integer :: ICLASS(maxconf,maxconf-1)
+ integer :: ICLASS(N,N-1)
+ INTEGER :: IA(N),IB(N),HVALS(LEV),IORDER(LEV),&
+ HEIGHT(LEV)
+ REAL(kind=4) :: CRIT(N),CRITVAL(LEV)
+ integer :: I,J,LOC,LEVEL,ICL,ILEV,NCL,K
+!
+! 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 (iout,450) (j,j=2,LEV)
+ 450 FORMAT(4X,' SEQ NOS',8(i2,'CL'),10000(i3,'CL'))
+ WRITE (iout,470) (' ---',j=2,LEV)
+ 470 FORMAT(4X,' -------',10000a4)
+ DO 500 I=1,N
+ WRITE (iout,600) I,(ICLASS(I,J),J=1,LEV-1)
+ 600 FORMAT(I11,8I4,10000i5)
+ 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 SUBROUTINE HCASS
+!-----------------------------------------------------------------------------
+!+++++++++++++++++++++++++++++++++++++++++++++++++C
+! C
+! 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)
+! include 'COMMON.IOUNITS'
+ integer :: LEV
+ CHARACTER(len=80) :: LINE
+ INTEGER :: IORDER(LEV),HEIGHT(LEV)
+ REAL(kind=4) :: CRITVAL(LEV)
+! INTEGER OUT(3*LEV,3*LEV)
+! INTEGER UP,ACROSS,BLANK
+ CHARACTER(len=1) :: OUT(3*LEV,3*LEV)
+ CHARACTER(len=1) :: UP,ACROSS,BLANK
+ DATA UP,ACROSS,BLANK /'|','-',' '/
+ integer :: I,I2,J,J2,K,I3,L,IC,IDUM
+!
+!
+ DO I=1,3*LEV
+ DO J=1,3*LEV
+ OUT(I,J)=BLANK
+ ENDDO
+ ENDDO
+!
+!
+ DO I=3,3*LEV,3
+ I2=I/3
+!
+ J2=3*LEV+1-3*HEIGHT(I2)
+ DO J=3*LEV,J2,-1
+ OUT(J,I)=UP
+ ENDDO
+!
+ DO K=I,3,-1
+ I3=INT((K+2)/3)
+ IF ( (3*LEV+1-HEIGHT(I3)*3).LT.J2) GOTO 100
+ OUT(J2,K)=ACROSS
+ ENDDO
+ 100 CONTINUE
+!
+ ENDDO
+!
+!
+ IC=3
+ DO I=1,3*LEV
+ IF (I.EQ.IC+1) THEN
+ IDUM=IC/3
+ IDUM=LEV-IDUM
+ DO L=1,LEV
+ IF (HEIGHT(L).EQ.IDUM) GOTO 190
+ ENDDO
+ 190 IDUM=L
+ WRITE(iout,200) CRITVAL(IDUM),(OUT(I,J),J=1,3*LEV)
+ IC=IC+3
+ ELSE
+ LINE = ' '
+ WRITE(iout,210) (OUT(I,J),J=1,3*LEV)
+ ENDIF
+ 200 FORMAT(1H ,8X,F12.2,4X,27000A1)
+ 210 FORMAT(1H ,24X,27000A1)
+ ENDDO
+ WRITE(iout,250)
+ WRITE(iout,220)(IORDER(J),J=1,LEV)
+ WRITE(iout,250)
+ 220 FORMAT(1H ,24X,9000I3)
+ WRITE(iout,230) LEV
+ 230 FORMAT(1H ,13X,'CRITERION CLUSTERS 1 TO ',i3)
+ WRITE(iout,240) LEV-1
+ 240 FORMAT(1H ,13X,'VALUES. (TOP ',i3,' LEVELS OF HIERARCHY).')
+ 250 FORMAT(/)
+!
+!
+ RETURN
+ END SUBROUTINE HCDEN
+!-----------------------------------------------------------------------------
+ end module hc_
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
projects / unres4.git / blobdiff