C @(#)clussbr.for	17.1.1.1 (ES0-DMD) 01/25/02 17:16:37
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION:
C      CLUS VERSION 1                 85.01.21
C      F. MURTAGH     ST-ECF, MUNICH
C
C.KEYWORDS:
C      APPLICATIONS, STATISTICS, CLUSTERING.
C
C.PURPOSE:
C      HIERARCHICAL CLUSTERING USING MINIMUM VARIANCE CRITERION.
C      N * M ARRAY PASSED TO PROGRAM IN "DATA",
C      SEQUENCE OF N-1 AGGLOMERATIONS OF ROW POINTS RETURNED IN
C      VECTORS "IA" AND "IB", WITH CRITERION VALUES IN "CRIT".
C
C.ALGORITHM:
C      NEAREST NEIGHBOUR CHAIN ALGORITHM IS USED (SEE BIBLIOGRA-
C      PHIC REFERENCES IN DOCUMENTATION), WHICH ALLOWS HIERARCHY
C      TO BE CONSTRUCTED WITHOUT THE NEED FOR A DISSIMILARITY
C      MATRIX.
C
C.INPUT/OUTPUT:
C      PARAMETERS:
C      N, M, DATA : DIMENSION OF, AND ARRAY, OF REALS;
C      IA, IB, CRIT : SEQ. OF AGGLOMERANDS & CRITERION VALUES,
C                     RETURNED (ONLY LOCNS 1 TO N-1 OF INTEREST);
C      MEMBR : CARDINALITIES OF CLASSES;
C      DISS  : NEAREST NEIGHBOUR DISSIMILARITIES;
C      ICHAIN: NEAREST NEIGHBOUR CHAIN;
C      FLAG  : INDICATES AGGLOMERABLE OBJECTS/CLUSTERS.
C      *** NOTE *** THAT STORAGE IS OVERLAID AS INDICATED BY
C                   COMMENTED EQUIVALENCE STATEMENT.
C                   ALSO NOTE THAT IJKLMN VARIABLE DEFINITION
C                   CONVENTION IS FOLLOWED.
C
C---------------------------------------------------------------
C
      SUBROUTINE CLUS(N,M,DATA,IA,IB,CRIT,
     X      MEMBR,DISS,ICHAIN,FLAG)
      REAL MEMBR(N),DATA(N,M),DISS(N),CRIT(N),D
      INTEGER ICHAIN(N),IA(N),IB(N),N,M,I,LEN,IDUM,I1,I2,NCL
      REAL INF
      LOGICAL FLAG(N)
      DATA INF/1.E+25/
C     EQUIVALENCE (ICHAIN(1),IA(1)),(DISS(1),CRIT(1))
C
C
      DO 150 I=1,N
      MEMBR(I)=1
      FLAG(I)=.TRUE.
  150 CONTINUE
      NCL=N
      I1=1
C
C-----START THE NN-CHAIN:--------------------------------------------
  200 LEN=N
      ICHAIN(LEN)=I1
      DISS(LEN)=INF
C-----DET. NN OF OBJECT I1:------------------------------------------
  300 FLAG(I1)=.FALSE.
C-----(TURN OFF FLAG SO THAT ZERO DISS. OF I1 WITH SELF NOT OBTAINED.)
      D=DISS(LEN)
      IF (LEN.LT.N) I2=ICHAIN(LEN+1)
C-----(FOR IDENTICAL DISS'S, ABOVE ENSURES THAT RNN WILL BE FOUND.)
      CALL DETNN(DATA,FLAG,MEMBR,N,M,I1,I2,D)
      FLAG(I1)=.TRUE.
C-----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
C-----IF LEN<N DISTINGUISH BETWEEN HAVING RNN & CONTINUING NN-CHAIN:--
  350 CONTINUE
      IF (I2.NE.ICHAIN(LEN+1)) GOTO 400
C-----(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
C-----(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
C-----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.
      IDUM=ICHAIN(LEN+1)
      IF (LEN.LT.N) FLAG(IDUM)=.TRUE.
C-----(REESTABLISH AGGLOMERABILITY OF OBJ.'S IN NN-CHAIN.)-------------
      GOTO 300
  550 CALL NEXT(FLAG,I1,N)
      GOTO 200
C
C
  600 CONTINUE
C     CALL OUTP(IA,IB,CRIT,N,M,DATA,dismax)
      RETURN
  700 WRITE(6,750)
  750 FORMAT(' ERROR IN NN-CHAIN ROUTINE - INSUFFICIENT CHAIN SPACE'/)
      STOP
      END
C
C
C
C
C
C-----------------------------------------------------------------------
C
C
C
C
C
      SUBROUTINE DETNN(DATA,FLAG,MEM,N,M,I1,I2,D)
C-----DETERMINE A NEAREST NEIGHBOUR; CALLED FROM SUBR. "CLUS".
      INTEGER N,M,I,J,I1,I2
      REAL DATA(N,M),MEM(N),D,DISS
      LOGICAL FLAG(N)
C
      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
C
      RETURN
      END
C
C
C
C
C
C-----------------------------------------------------------------------
C
C
C
C
C
      SUBROUTINE AGGLOM(I1,I2,D,DATA,MEM,FLAG,IA,IB,CRIT,NCL,N,M)
C-----CARRY OUT AN AGGLOMERATION; CALLED FROM SUBR. "CLUS".
      REAL MEM(N),DATA(N,M),CRIT(N),D
      INTEGER IA(N),IB(N),N,M,NCL,I1,I2,I,J,NAGGL
      LOGICAL FLAG(N)
      INTEGER O1,O2,LB,UB,J1
C
C
      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) )
     X            / (MEM(O1)+MEM(O2))
      DATA(O2,J)=DATA(O1,J)
  100 CONTINUE
      NAGGL=N-NCL
      MEM(O1)=MEM(O1)+MEM(O2)
      FLAG(O2)=.FALSE.
C-----KEEP SORTED LIST OF CRIT. VALUES: FIND 1ST WHERE NEW CRIT. FITS:
      I=NAGGL-1
  120 IF (D.GE.CRIT(I)) GOTO 140
      I=I-1
      IF (I.GE.1) GOTO 120
C-----(ARRIVING HERE MUST MEAN THAT D > ALL CRIT. VALUES FOUND SO FAR:)
      I=0
  140 CONTINUE
C-----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
C
C
      RETURN
      END
C
C
C
C
C
C-----------------------------------------------------------------------
C
C
C
C
C
      SUBROUTINE NEXT(FLAG,I1,N)
C-----DETERMINE NEXT AGGLOMERABLE OBJECT/CLUSTER; CALLED FROM "CLUS".
      INTEGER N,I1,NXT,I
      LOGICAL FLAG(N)
C
C
      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
C
C
      STOP
C
C
  500 I1=I
C
C
      RETURN
      END
C
C
C
C
C
C-----------------------------------------------------------------------
C
C
C
C
C
      SUBROUTINE OUTP(IA,IB,CRIT,N,M,DATA)
C-----OUTPUT DATA: THIS ROUTINE IS CURRENTLY DEAD.
      INTEGER N,M,IA(N),IB(N),NMIN1,JJJ,K
      REAL CRIT(N),DATA(N,M)
C
C
      NMIN1=N-1
      DO 100 JJJ=1,NMIN1
      WRITE(6,50) JJJ,IA(JJJ),IB(JJJ),CRIT(JJJ)
   50 FORMAT(1H ,I6,I6,I6,F15.4)
  100 CONTINUE
C
C
      WRITE(6,70)
   70 FORMAT(////,' CLUSTER CENTRE COORDINATES FOLLOW.',///)
      DO 300 JJJ=1,N
      WRITE(6,80) JJJ,(DATA(JJJ,K),K=1,M)
   80 FORMAT(I4,13F6.3)
  300 CONTINUE
C
      RETURN
      END