!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),DISS(18) ! REAL CRITVAL(9) ! INTEGER IA(18),IB(18),ICHAIN(18),ICLASS(18,9),HVALS(9) ! INTEGER IORDER(9),HEIGHT(9) ! LOGICAL FLAG(18) !Integer, parameter :: N=36, M=7, NN2=N*(N-1)/2, Lev=5 !REAL DATA(N,M),CRIT(N),MEMBR(N),DISS(N) !REAL CRITVAL(Lev) !INTEGER IA(N),IB(N),ICHAIN(N),ICLASS(N,Lev),HVALS(Lev) !INTEGER IORDER(Lev),HEIGHT(Lev) !LOGICAL FLAG(N) ! IN ABOVE, 18=N, 16=M, 9=LEV, 153=N(N-1)/2. ! Use NetCDF Integer :: N, M, NN2, Lev Real, Allocatable :: DATA(:,:),CRIT(:),MEMBR(:),DISS(:) Real, Allocatable :: CRITVAL(:) Integer, Allocatable :: IA(:),IB(:),ICHAIN(:),ICLASS(:,:),HVALS(:) Integer, Allocatable :: IORDER(:),HEIGHT(:) Logical, Allocatable :: FLAG(:) Integer :: Nst, Ntime Real, Allocatable :: StnID(:), X(:,:), Lon(:), Lat(:), Reg(:), xtime(:) ! ! 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) Call Read_PDSI !N=Nst; M=Ntime; Lev=5 N=Nst; M=2; Lev=7 Allocate (DATA(N,M),CRIT(N),MEMBR(N),DISS(N)) Allocate (CRITVAL(Lev)) Allocate (IA(N),IB(N),ICHAIN(N),ICLASS(N,Lev),HVALS(Lev)) Allocate (IORDER(Lev),HEIGHT(Lev),FLAG(N)) Data(:,:)=X(:,:) !Data(1:N,1)=Lon(1:N); Data(1:N,2)=Lat(1:N) CALL HCON2(N,M,DATA,IA,IB,CRIT,MEMBR,DISS,ICHAIN,FLAG) CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT,Reg) !CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL) !-------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------- Subroutine Read_PDSI Character (Len=200) :: Fname Fname="/home/tanpv/SEACLID_New/RES/PDSI_Res/"// & "PDSI_Stn_No_Rain_Center.nc" !Call rd_ncdf_2d (Fname,StnID,Xtime,Nst,Ntime,X,"ori_pdsi") Call rd_ncdf_2d (Fname,StnID,Xtime,Nst,Ntime,X,"Pre") Deallocate (StnID) Call rd_ncdf_1d (Fname, StnID, Nst, Lon, "lon") Deallocate (StnID) Call rd_ncdf_1d (Fname, StnID, Nst, Lat, "lat") Deallocate (StnID) Call rd_ncdf_1d (Fname, StnID, Nst, Reg, "reg") print*, Reg End Subroutine Read_PDSI !-------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" !-------------------------------------------------------------------- END !--------------------------------------------------------------- ! HIERARCHICAL CLUSTERING using Minimum Variance Criterion, C ! using the O(N**2) time Nearest Neighbour Chain algorithm. C ! C ! Parameters: C ! C ! DATA(N,M) : input data, C ! IA(N), IB(N), CRIT(N) : sequence of agglomerands and C ! values returned (only locations 1 to N-1 C ! are of interest), C ! MEMBR(N), DISS(N), ICHAIN(N) : used in the routines to store C ! cluster cardinalities, nearest neighbour C ! dissimilarities, and the NN-chain. C ! FLAG(N) : (boolean) used to indicate agglomerable C ! objects and clusters. C ! C ! Reference: Murtagh, Multidimensional Clustering Algorithms, C ! Physica-Verlag, 1985. C ! C ! F. Murtagh, ESA/ESO/STECF, Garching, February 1986. C ! C ! HISTORY C ! C ! Bounds bug fix, Oct. 1990, F. Murtagh. C ! On line 98, alter stmt. "IDUM=ICHAIN(LEN+1)" to C ! "IF (LEN.LT.N) IDUM=ICHAIN(LEN+1)" C ! The flag in the following stmt. took care of this eventuality C ! and the VAX/VMS compiler, without the /CHECK=(BOUNDS) option C ! did not care about the fact that ICHAIN went out of bounds. C !----------------------------------------------------------------C SUBROUTINE HCON2(N,M,DATA,IA,IB,CRIT,MEMBR,DISS,ICHAIN,FLAG) REAL MEMBR(N), DATA(N,M), DISS(N), CRIT(N) INTEGER ICHAIN(N), IA(N), IB(N) REAL INF LOGICAL FLAG(N) DATA INF/1.E+25/ ! EQUIVALENCE (ICHAIN(1),IA(1)),(DISS(1),CRIT(1)) DO 150 I=1,N MEMBR(I)=1 FLAG(I)=.TRUE. 150 CONTINUE NCL=N I1=1 ! Start the NN-chain: 200 LEN=N ICHAIN(LEN)=I1 DISS(LEN)=INF ! Determine NN of object I1 300 FLAG(I1)=.FALSE. ! Turn off FLAG so that 0 diss. of I1 with self not obtained. D=DISS(LEN) IF (LEN.LT.N) I2=ICHAIN(LEN+1) ! For identical diss.'s, above ensures that RNN will be found. CALL DETNN(DATA,FLAG,MEMBR,N,M,I1,I2,D) FLAG(I1)=.TRUE. ! If LEN = 1 place obj. I2 as second obj. in NN-chain. IF (LEN.LT.N) GOTO 350 LEN=LEN-1 IF (LEN.LT.N-NCL) GOTO 700 ICHAIN(LEN)=I2 DISS(LEN)=D GOTO 500 ! If LEN < N distinguish between having RNN & continuing NN-chain. 350 CONTINUE IF (I2.NE.ICHAIN(LEN+1)) GOTO 400 ! Have RNN. NCL=NCL-1 CALL AGGLOM(I1,I2,D,DATA,MEMBR,FLAG,IA,IB,CRIT,NCL,N,M) LEN=LEN+2 GOTO 500 400 CONTINUE ! Grow extra link on NN-chain. IDUM=ICHAIN(LEN+1) FLAG(IDUM)=.FALSE. LEN=LEN-1 IF (LEN.LE.N-NCL) GOTO 700 ICHAIN(LEN)=I2 DISS(LEN)=D GOTO 500 ! Select obj. for continuing to grow (or restarting) NN-chain. 500 CONTINUE IF (NCL.EQ.1) GOTO 600 IF (LEN.EQ.N+1) GOTO 550 I1=ICHAIN(LEN) FLAG(I1)=.TRUE. IF (LEN.LT.N) IDUM=ICHAIN(LEN+1) IF (LEN.LT.N) FLAG(IDUM)=.TRUE. ! Reestablish agglomerability of objects in NN-chain. GOTO 300 550 CALL NEXT(FLAG,I1,N) GOTO 200 600 CONTINUE RETURN 700 WRITE(6,750) 750 FORMAT(' ERROR IN NN-CHAIN ROUTINE - INSUFFICIENT CHAIN SPACE'/) STOP END !-------------------------------------------------------- SUBROUTINE DETNN(DATA,FLAG,MEM,N,M,I1,I2,D) ! Determine a nearest neighbour. ! REAL DATA(N,M),MEM(N) LOGICAL FLAG(N) DO 200 I=1,N IF (.NOT.FLAG(I)) GOTO 200 DISS=0. DO 100 J=1,M 100 DISS=DISS+(DATA(I1,J)-DATA(I,J))*(DATA(I1,J)-DATA(I,J)) DISS=DISS*MEM(I)*MEM(I1)/(MEM(I1)+MEM(I)) IF (DISS.GE.D) GOTO 200 D=DISS I2=I 200 CONTINUE RETURN END !------------------------------------------------------- SUBROUTINE AGGLOM(I1,I2,D,DATA,MEM,FLAG,IA,IB,CRIT,NCL,N,M) ! HISTORY ! Bounds bug fix, Oct. 1990, F. Murtagh. ! Insert line 25, "IF (I.EQ.0) GOTO 140". Without /CHECK=(BOUNDS) option ! the VAX-VMS compiler did not complain. ! ! Carry out an agglomeration. ! REAL MEM(N),DATA(N,M),CRIT(N) INTEGER IA(N),IB(N) LOGICAL FLAG(N) INTEGER O1,O2,LB,UB O1=MIN0(I1,I2) O2=MAX0(I1,I2) DO 100 J=1,M DATA(O1,J)=( MEM(O1)*DATA(O1,J)+MEM(O2)*DATA(O2,J) ) / (MEM(O1)+MEM(O2)) DATA(O2,J)=DATA(O1,J) 100 CONTINUE NAGGL=N-NCL MEM(O1)=MEM(O1)+MEM(O2) FLAG(O2)=.FALSE. ! Keep sorted list of criterion values: find 1st where new crit. fits. I=NAGGL-1 IF (I.EQ.0) GOTO 140 120 IF (D.GE.CRIT(I)) GOTO 140 I=I-1 IF (I.GE.1) GOTO 120 ! Arriving here must mean that D > all crit. values found so far. I=0 140 CONTINUE ! Now, shift rightwards from I+1 to AGGL-1 to make room for new crit. LB=I+1 UB=NAGGL-1 IF (LB.GT.UB) GOTO 180 J=UB 160 J1=J+1 IA(J1)=IA(J) IB(J1)=IB(J) CRIT(J1)=CRIT(J) J=J-1 IF (J.GE.LB) GOTO 160 180 CONTINUE IA(LB)=O1 IB(LB)=O2 CRIT(LB)=D RETURN END !------------------------------------------- SUBROUTINE NEXT(FLAG,I1,N) ! Determine next agglomerable object/cluster. LOGICAL FLAG(N) NXT=I1+1 IF (NXT.GT.N) GOTO 150 DO 100 I=NXT,N IF (FLAG(I)) GOTO 500 100 CONTINUE 150 DO 200 I=1,I1 IF (FLAG(I)) GOTO 500 200 CONTINUE STOP 500 I1=I 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,Reg) INTEGER IA(N),IB(N),ICLASS(N,LEV),HVALS(LEV),IORDER(LEV),HEIGHT(LEV) REAL CRIT(N),CRITVAL(LEV),Reg(N) ! 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,Int(Reg(I)),(ICLASS(I,J),J=1,Lev) !WRITE (6,600) I,(ICLASS(I,J),J=1,Lev) 600 FORMAT(I11,15I10) 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 !--------------------------------------------------------------- !---------------------------------------------------------------