*-----------------------------------------------------------------------
* cluster.f - data clustering for music analysis
* by Andy Allinger, 2011-2017, released to the public domain
* This program may be used by any person for any purpose.
*-----------------------------------------------------------------------

*-----------------------------------------------------------------------
* SMD - Subtract Mixed Data
*
*__Name______Type______In/Out_________Description_______________________
*  A         Real      In             Value
*  B         Real      In             Value
*  T         Real      In             Data type descriptor
*  SMD       Real      Out            Absolute difference
*-----------------------------------------------------------------------
      FUNCTION SMD (A, B, T)       
       IMPLICIT NONE
       REAL A, B, T, SMD
       
       REAL LIT, RMACHINE
       LOGICAL FIRST
       SAVE LIT, FIRST
       DATA FIRST / .TRUE. /
       IF (FIRST) THEN
         LIT = RMACHINE(2)
         FIRST = .FALSE.
       END IF
       
       IF (ABS(T) < LIT) THEN                   ! common math
         SMD = ABS(A - B)
       ELSE IF (T > 0.) THEN                        ! modulo arithmetic
         SMD = MOD( ABS(A - B), T )
         SMD = MIN(SMD, T - SMD)        ! choose the shorter path
       ELSE IF (T > -1.5) THEN             ! categorical difference
         IF (ABS(A - B) > 0.5) THEN                 
           SMD = 1.
         ELSE
           SMD = 0.
         END IF
       ELSE                             ! contra-categorical
         IF (ABS(A - B) > 0.5) THEN
           SMD = 0.
         ELSE
           SMD = 1.
         END IF
       END IF
       RETURN
      END ! of smd


*-----------------------------------------------------------------------
* DIST1 - the Manhattan distance between two points 
*         D = SUM (W[i] @ |(A[i]-B[i]) % T[i]|)  
*
*__Name_____Type_____In/Out___Description_______________________________
*  A[P]     Real     In       Column of a data matrix
*  B[P]     Real     In       Column of a data matrix
*  T[P]     Real     In       Modulo arithmetic divisors
*  W[P]     Real     In       Importance weights
*  P        Integer  In       Dimension of space
*  DIST1    Real     Out      Manhattan distance
*-----------------------------------------------------------------------
      FUNCTION DIST1 (A, B, T, W, p)
       IMPLICIT NONE
       INTEGER p
       REAL A, B, T, W, DIST1
       DIMENSION A(P), B(P), T(P), W(P)
                   
       INTEGER J
       REAL DIF, SMD

       DIST1 = 0.
       DO J = 1, P
         DIF = SMD(A(J), B(J), T(J))
         DIST1 = DIST1 + W(J) * DIF
       END DO
       RETURN
      END  ! of DIST1
      
      
*-----------------------------------------------------------------------
* DIST2 - the squared Euclidean distance between two points 
*         D = SUM (W[i] @ (A[i]-B[i]) % T[i]) ^2 
*__Name_____Type_____In/Out___Description_______________________________
*  A[P]     Real     In       Column of a data matrix
*  B[P]     Real     In       Column of a data matrix
*  T[P]     Real     In       Modulo arithmetic divisors
*  W[P]     Real     In       Importance weights
*  P        Integer  In       Dimension of space
*  DIST2    Real     Out      Squared Euclidean distance
*-----------------------------------------------------------------------
      FUNCTION DIST2(A, B, T, W, P)
       IMPLICIT NONE
       INTEGER P
       REAL A, B, T, W, DIST2
       DIMENSION A(P), B(P), T(P), W(P)

       INTEGER J
       REAL DIF, SMD

       DIST2 = 0.
       DO J = 1, P
         DIF = SMD(A(J), B(J), T(J))
         DIST2 = DIST2 + (W(J) * DIF) **2
       END DO
       RETURN
      END     ! of DIST2


*=======================================================================
* hclust - Hyperplane partitioning
* Divisive clustering based on partial sorting, a reduced form of
* Anthony Hoare's quick sort algorithm.
*
*_____Name_____________Type________I/O___Description____________________
*     X[PLIM,N]        REAL        In    data array
*     PLIM             INTEGER     In    row dimension of X
*     P                INTEGER     In    dimension of data
*     N                INTEGER     In    number of data points
*     T[P]             REAL        In    modulo period for distance
*     W[P]             REAL        In    importance weights
*     K                INTEGER     In    desired number of clusters
*     Z[N]             INTEGER     Out   what cluster each point is in
*     POP[K]           REAL        Out   cluster populations
*     IND[N]           INTEGER     Out   index array that sorts X into
*                                         output order from input order
*     BOUND[K]         INTEGER     None  work array
*     DFUNC            FUNCTION    In    function name of return type REAL
*                                         dfunc(A, B, T, W, p)
*                                         INTEGER p
*                                         REAL A[p], B[p], T[p], W[p]
*     FAULT            INTEGER     Out   exit code
*=======================================================================
      SUBROUTINE HCLUST
     &      (X,PLIM,P,N,T,W,K,Z,POP,IND,BOUND,DFUNC,FAULT)

       IMPLICIT NONE
       INTEGER PLIM, P, N, K, Z, IND, BOUND, FAULT
       REAL X, T, W, POP, DFUNC
       DIMENSION X(PLIM,N), T(P), W(P),
     &    Z(N), POP(K), IND(N), BOUND(K)

*          Local variables
       INTEGER HI, LO, IHI, ILO, i, L, LL, HAT
       INTEGER  BIGSIZE, BIGPART, PARTSIZE, ISWAP, RL, RM, RH
       REAL DHL, DHM, DML

*         initialize
       FAULT = 0
       DO i = 1, N
         IND(i) = i
       END DO
       DO L = 1, K
         BOUND(L) = N
       END DO
       POP(1) = FLOAT(N)
       
*          Trivial cases
       IF (K .EQ. 1) THEN
         DO i = 1, N
           Z(i) = 1
         END DO
         POP(1) = FLOAT(N)
         RETURN
       ELSE IF (K .EQ. N) THEN
         DO i = 1, N
           Z(i) = i
         END DO
         DO L = 1, K
           POP(L) = 1.
         END DO
         RETURN
       END IF

*=======================================================================
*       Partition a section of the array such that elements closer to the first
*       element are moved to the front, and elements closer to the last element
*       are moved to the rear.  Continue until enough clusters have been made.
*=======================================================================
       DO L = 1, K-1

         PARTSIZE = NINT(POP(1))          ! find largest partition
         BIGSIZE = PARTSIZE
         BIGPART = 1
         DO LL = 2, L
           PARTSIZE = NINT(POP(LL))
           IF (PARTSIZE > BIGSIZE) THEN
             BIGSIZE = PARTSIZE
             BIGPART = LL
           END IF
         END DO

         HI = BOUND(BIGPART)               ! find limits of partition
         LO = HI - BIGSIZE + 1
         IHI = HI - 1
         ILO = LO + 1
         IF (LO < 1) PRINT *, '!!lo=', lo, ' bigpart=', bigpart, bigsize
         IF (BIGSIZE .EQ. 2) GO TO 70

*          swap random pivot values to ends of partition (introduce chance)
         RH = HAT(HI-LO+1) -1 + LO
  20     RL = HAT(HI-LO+1) -1 + LO
         IF (RL .EQ. RH) GO TO 20

*           a third pivot value in case first two are near each other
  30     RM = HAT(HI-LO+1) -1 + LO
         IF (RM .EQ. RH .OR. RM .EQ. RL) GO TO 30
         DHL = DFUNC(X(1,IND(RH)), X(1,IND(RL)), T, W, P)
         DHM = DFUNC(X(1,IND(RH)), X(1,IND(RM)), T, W, P)
         DML = DFUNC(X(1,IND(RM)), X(1,IND(RL)), T, W, P)
         IF (DHM > DHL .OR. DML > DHL) THEN
           IF (DHM > DML) THEN
             RL = RM
           ELSE
             RH = RM
           END IF
         END IF

         ISWAP = IND(RH)
         IND(RH) = IND(HI)
         IND(HI) = ISWAP
         ISWAP = IND(RL)
         IND(RL) = IND(LO)
         IND(LO) = ISWAP

*=======================================================================
*            sort
*=======================================================================
  50     IF (ILO .LE. IHI) THEN       ! compare to lo pivot
           IF (DFUNC(X(1,IND(ILO)), X(1,IND(LO)), T, W, P) .LE.
     &    DFUNC(X(1,IND(ILO)), X(1,IND(HI)), T, W, P)) THEN
             ILO = ILO + 1
             GO TO 50
           END IF
         END IF

  60     IF (IHI .GE. ILO) THEN        ! compare to hi pivot
           IF (DFUNC(X(1,IND(IHI)), X(1,IND(HI)), T, W, P) .LE.
     &    DFUNC(X(1,IND(IHI)), X(1,IND(LO)), T, W, P)) THEN
             IHI = IHI - 1
             GO TO 60
           END IF
         END IF

         IF (ILO < IHI) THEN
           ISWAP = IND(ILO)
           IND(ILO) = IND(IHI)
           IND(IHI) = ISWAP
           ILO = ILO + 1
           IHI = IHI - 1
           IF (ILO .LE. IHI) GO TO 50
         END IF

*=======================================================================
  70   CONTINUE      ! Add this partition boundary to the bounds array
*=======================================================================
*             equal split if all distances were the same
         IF (ILO .EQ. HI) THEN
           IHI = (LO + HI) / 2
           ILO = IHI + 1
         END IF
         POP(BIGPART) = FLOAT(HI - IHI)        ! reduce population
         POP(L+1) = FLOAT(ILO - LO)            ! population in new partition
         BOUND(L+1) = IHI               ! append
       END DO     ! next partition

*=======================================================================
*  Recover memberships and populations from the index array
*=======================================================================
       CALL ISORT (BOUND, K)
       L = 1
       DO i = 1, N
         IF (i > BOUND(L)) L = L + 1
         Z(IND(i)) = L
       END DO

       i = 0
       DO L = 1, K
         POP(L) = FLOAT(BOUND(L) - i)
         i = BOUND(L)
       END DO
       RETURN
      END  ! of HCLUST



*-----------------------------------------------------------------------
* Passive sort of unsigned integers by an unstable radix algorithm
*
*___Name_____Type______In/Out___Description_____________________________
*   IX(N)    Integer   In       Input array
*   IND(N)   Integer   Out      Permutation vector
*   N        Integer   In       Length of array
*   M        Integer   In       Largest value in array
*-----------------------------------------------------------------------
      SUBROUTINE USORTU (IX, IND, N, M)
       IMPLICIT NONE
       INTEGER IX, IND, N, M
       DIMENSION IX(N), IND(N)

       INTEGER B, I, J, K, L, IL(0:31), IG(0:31), LG(0:31)  ! local variables
       INTEGER BITSZ                                        ! external function

*         Begin
       DO I = 1, N                ! initialize
         IND(I) = I
       END DO
       IF (N < 2) RETURN          ! validate

       J = BITSZ() - 1             ! greatest bit position
       B = J                       ! avoid compiler warning
       DO I = J, 0, -1
         B = I
         IF (BTEST(M, I)) GO TO 10
       END DO

  10   L = B                 ! most significant bit
       IL(L) = 1
       IG(L) = N
       LG(L) = -1

*            Sort
  20   I = IL(L)
       J = IG(L)
       IF (I .GE. J) GO TO 50    ! if empty or single, back up

*           Find a set bit at the front of the partition
  30   IF (I .LE. J) THEN
         IF (.NOT. BTEST(IX(IND(I)), L)) THEN
           I = I + 1
           GO TO 30
         END IF
       END IF

*           Find a clear bit at the rear of the partition
  40   IF (I .LE. J) THEN
         IF (BTEST(IX(IND(J)), L)) THEN
           J = J - 1
           GO TO 40
         END IF
       END IF

*             Swap these elements
       IF (I < J) THEN
         K = IND(I)
         IND(I) = IND(J)
         IND(J) = K
         I = I + 1
         J = J - 1
         GO TO 30
       END IF

*              Traverse tree
       IF (L .EQ. 0) GO TO 50  ! At the least significant bit, go back up
       IL(L-1) = IL(L)
       IG(L-1) = J
       LG(L) = -1
       L = L - 1                     ! Descend
       GO TO 20

  50   L = L + 1                     ! Back up
       IF (L > B) RETURN
       IF (LG(L) > 0) GO TO 50  ! Both subpartitions sorted; go back up
       IL(L-1) = IG(L-1) + 1         ! Do the greater sub-region
       IG(L-1) = IG(L)
       LG(L) = +1                    ! Mark as greater
       L = L - 1                     ! Descend
       GO TO 20
      END   ! of USORTU



*-----------------------------------------------------------------------
*    Let each center be the average of its points
*
*___Name____________Type______In/Out____Description_____________________
*   X(PLIM,N)       Real      In        Data points
*   PLIM            Integer   In        Row dimension of X
*   P               Integer   In        Variables each datum
*   N               Integer   In        # of data points
*   T(P)            Real      In        Periods
*   CAT(P)          Integer   In        # of categories
*   CATTOT          Integer   In        Sum of CAT()
*   C(PLIM,K)       Real      Out       Centers
*   K               Integer   In        Number of clusters
*   TAL(CATTOT,K)   Real      Out       Frequencies of observers
*   Z(N)            Integer   In        Cluster membership
*   POP(K)          Real      Out       Cluster population
*   ROBUST          Logical   In        Median instead of mean
*   WANTM           Logical   In        Compute modes of categoricals
*   LIVEC(K)        Integer   In        Skip clusters with livec .LE. 0
*   IPERM(N)        Integer   Neither   Workspace
*   XW(2@N)         Real      Neither   Workspace
*-----------------------------------------------------------------------
      SUBROUTINE CENTER (X, PLIM, P, N, T, CAT, CATTOT,
     &  C, K, TAL, Z, POP, ROBUST, WANTM, LIVEC, IPERM, XW)

       IMPLICIT NONE
       INTEGER PLIM, P, N, CAT, CATTOT, K, Z, LIVEC, IPERM
       REAL X, T, C, TAL, POP, XW
       LOGICAL ROBUST, WANTM
       DIMENSION X(PLIM,N), T(P), CAT(P), C(PLIM,K),
     &    TAL(CATTOT,K), Z(N), POP(K), LIVEC(K), IPERM(N), XW(2*N)

*        Local variables
       INTEGER H, HCUM, I, I1, I5, I9, J, L
       REAL LIT, DEV
       LOGICAL FIRST

*         External functions
       INTEGER HAT
       REAL RMACHINE

       SAVE LIT, FIRST
       DATA FIRST / .TRUE. /

*-----------------------------------------------------------------------
*           begin
*-----------------------------------------------------------------------
       IF (FIRST) THEN
         LIT = RMACHINE(2)
         FIRST = .FALSE.
       END IF

*       Count cluster populations.
       DO L = 1, K
         POP(L) = 0.
       END DO
       DO I = 1, N
         L = Z(I)
         POP(L) = POP(L) + 1.
       END DO

*       Sort implicitly by cluster ID
       CALL USORTU (Z, IPERM, N, K)

*        averages
       DO 30 J = 1, P
         IF (CAT(J) > 0 .AND. .NOT. WANTM) GO TO 30    ! don't need mode
         I9 = 0
         DO 20 L = 1, K
           I1 = I9 + 1
           I9 = I9 + NINT(POP(L))
           IF (LIVEC(L) .LE. 0 .OR. POP(L) < 0.5) GO TO 20   ! skip dead/empty
           H = 0
           DO I5 = I1, I9
             I = IPERM(I5)
             H = H + 1
             XW(H) = X(J,I)
           END DO
           IF (CAT(J) > 0) THEN
             CALL RMODE (XW, H, CAT(J), TAL, I)
             I = HAT(I)
             C(J,L) = TAL(I,1)         ! OK because i .LE. cattot always
           ELSE IF (T(J) > LIT) THEN        ! modulo arithmetic
             IF (ROBUST) THEN
               CALL CYMED (XW, H, T(J), XW(H+1), I, DEV)
             ELSE
               CALL CYMEAN (XW, H, T(J), XW(H+1), I, DEV)
             END IF
             I = HAT(I)
             C(J,L) = XW(H+I)
           ELSE IF (ROBUST) THEN           ! scalar data
             CALL MEDIAN (XW, H, C(J,L))
           ELSE                             ! arithmetic mean
             CALL MEAN (XW, H, C(J,L))
           END IF   ! robust ?
  20     CONTINUE
  30   CONTINUE
       IF (WANTM) RETURN     ! don't need to tally in this case

*          categorical and contra-categorical frequencies
       DO L = 1, K                ! zero
         DO H = 1, CATTOT
           TAL(H,L) = 0.
         END DO
       END DO

       HCUM = 0
       DO 50 J = 1, P
         IF (CAT(J) .EQ. 0) GO TO 50     ! skip others
         DO I = 1, N                              ! count
           H = HCUM + NINT(X(J,I))
           L = Z(I)
           TAL(H,L) = TAL(H,L) + 1.
         END DO
         HCUM = HCUM + CAT(J)
  50   END DO  ! next j
       RETURN
      END   ! of center


*-----------------------------------------------------------------------
* Let each center be one of its points
*
*___Name____________Type______In/Out____Description_____________________
*   X(PLIM,N)       Real      In        Data points
*   PLIM            Integer   In        Row dimension of X
*   P               Integer   In        Variables each datum
*   N               Integer   In        # of data points
*   C(PLIM,K)       Real      Out       Centers
*   K               Integer   In        Number of clusters
*   Z(N)            Integer   In        Cluster membership
*   POP(K)          Real      Out       Cluster population
*   IPERM(N)        Integer   Neither   Workspace
*-----------------------------------------------------------------------
      SUBROUTINE XOVER (X, PLIM, P, N, C, K, Z, POP, IPERM)
       IMPLICIT NONE
       INTEGER PLIM, P, N, K, Z, IPERM
       REAL X, C, POP
       DIMENSION X(PLIM,N), C(PLIM,K), Z(N), POP(K), IPERM(N)

*        Local variables
       INTEGER I, I5, I9, J, L, NPL

*         External functions
       INTEGER HAT

*-----------------------------------------------------------------------
*           begin
*-----------------------------------------------------------------------
*       Count cluster populations.
       DO L = 1, K
         POP(L) = 0.
       END DO
       DO I = 1, N
         L = Z(I)
         POP(L) = POP(L) + 1.
       END DO

!        PRINT *, 'xover: input k is: ', k
!        K = 0
!        DO I = 1, N
!          K = MAX(K, Z(I))
!        END DO
!        PRINT *, 'xover: maximum Z(i) is: ', K

*       Sort implicitly by cluster ID
       CALL USORTU (Z, IPERM, N, K)
       
!       PRINT *, 'xover: iperm is now: ', IPERM

*       Choose a random member
       I9 = 0
       DO 20 L = 1, K
         NPL = NINT(POP(L))
         I5 = I9 + HAT(NPL)
         I9 = I9 + NPL
         IF (NPL < 1) GO TO 20
         I = IPERM(I5)
!         PRINT *, 'i is ', i, ' and i5 is ', i5
         DO J = 1, P
           C(J,L) = X(J,I)
         END DO
  20   CONTINUE
       RETURN
      END  ! of XOVER


*-----------------------------------------------------------------------
* recs - re-estimate by cubic smoothing
*
*___Name________Type______I/O_______Description_________________________
*   X(PLIM,N)   Real      Both      Data points in R^p
*   PLIM        Integer   In        Row dimension of X
*   P           Integer   In        Measurements per point
*   N           Integer   In        Number of points
*   T(P)        Real      In        Period of modulo arithmetic
*   M           Integer   In        Number of sources
*   C(PLIM,K)   Real      In        Cluster centroids
*   K           Integer   In        # clusters
*   Z(N)        Integer   In        What cluster each point is in
*   IPERM(N)    Integer   Neither   Permutation vector
*   IDX(N)      Integer   Neither   File ID, each point
*   RW(N)       Real      In        W values for CFIT
*   RX(N)       Real      Work      X values for CFIT
*   RY(N)       Real      Neither   Y values for CFIT
*-----------------------------------------------------------------------
      SUBROUTINE RECS (X,PLIM,P,N,T,M,C,K,Z,IPERM,IDX,RW,RX,RY)
       IMPLICIT NONE
       INTEGER PLIM, P, N, M, K, Z, IPERM, IDX
       REAL X, T, C, RW, RX, RY
       DIMENSION X(PLIM,N), T(P), C(PLIM,K), Z(N), 
     $  IPERM(N), IDX(N), RW(N), RX(N), RY(N)

*          Constants
       REAL SMOOTH
       PARAMETER (SMOOTH = 0.05)  ! amount to adjust towards regression value

*            Local variables
      INTEGER
     $        H, HOLD,          ! file ID's
     $        I, I0, I2, I5,    ! count points
     $        IERR,             ! error code
     $        J,                ! count dimensions
     $        L,                ! cluster index
     $        O                 ! # objects in input file
     
       REAL A(0:3),            ! cubic coefficients
     $      R                  ! temp scalar

*             External functions
!!       REAL CBRT                    ! cube root
       LOGICAL ALMEQ                ! the almost-equal function

*-----------------------------------------------------------------------
*        Begin.
*-----------------------------------------------------------------------
       IERR = 0                    ! clear error codes
       DO I = 1, N                 ! source file (refer to BLEND in mixmidi.f)
         IDX(I) = NINT(X(3,I))
       END DO
       CALL USORTU (IDX, IPERM, N, M)  ! Sort implicitly
       
       DO 10 J = 1, P              ! for each dimension
         IF (.NOT. ALMEQ(T(J), 0.0)) GO TO 10
         DO I = 1, N                ! fill work arrays
           I5 = IPERM(I)
           RX(I) = X(J,I5)
           L = Z(I5)
           RY(I) = C(J,L)
         END DO

         H = IDX(IPERM(1))           ! which source file
         I0 = 1
         O = 0
         DO I = 1, N                 ! for each object
           I5 = IPERM(I)
           HOLD = H
           H = IDX(I5)
           IF (I .EQ. N) O = O + 1
           IF (H .NE. HOLD .OR. I .EQ. N) THEN      ! compute cubic
           
!!!!!!!!!!!!!!!!!!!! debug
!             PRINT *, 'data sent to cfit:'
!             DO I2 = I0, I0+O-1
!               PRINT *, RW(I2),  RX(I2), RY(I2)
!             END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           
             CALL CFIT (RW(I0), RX(I0), RY(I0), O, A, IERR)
             IF (IERR .NE. 0) PRINT *, 'recs: cfit returns #', IERR
             DO I2 = I0, I0+O-1
               R = RX(I2)
               RY(I2) = ((A(3) * R + A(2)) * R + A(1)) * R + A(0)
             END DO  ! next i
             I0 = I
             O = 0
           END IF
           O = O + 1
         END DO  ! next i
         
*             Apply the estimate
   8     FORMAT ('file ', I3, ' changing X from: ', F14.6, $)
   9     FORMAT (' to: ', F14.6)
         DO I = 1, N
           I5 = IPERM(I)
!           PRINT 8, IDX(I5), X(J,I5)
            X(J,I5) = X(J,I5) + SMOOTH * (RY(I) - X(J,I5))
!           PRINT 9, X(J,I5)
         END DO         
  10   CONTINUE  ! next j
       RETURN
      END  ! of RECS


************************************************************************
* DSM - a dissimilarity function
*
*___Variable______Type______I/O_______Description_______________________
*   X(P)          Real      In        Data points
*   P             Integer   In        Dimension of the data
*   T(P)          Real      In        Period of modulo arithmetic
*   W(P)          Real      In        Importance weights
*   CAT(P)        Integer   In        # categories each variable
*   CATTOT        Integer   In        Total categories
*   C(P)          Real      Both      Center points of clusters
*   TAL(CATTOT)   Real      Both      Observer frequencies
*   POP           Real      In        # of points in each cluster
*   LIT           Real      In        Negliglible value
*   SELF          Real      In        Object is in this cluster?
*   ROBUST        Logical   In        Use medians
*   BEST          Real      In        Bail out if this value is exceeded
*   DSM           Real      Out       Dissimilarity
************************************************************************
      FUNCTION DSM (X,P,T,W,CAT,CATTOT,C,TAL,POP,LIT,SELF,ROBUST,BEST)
       IMPLICIT NONE
       INTEGER P, CAT, CATTOT
       REAL X, T, W, C, TAL, POP, LIT, SELF, BEST, DSM
       LOGICAL ROBUST
       DIMENSION X(P), T(P), W(P), CAT(P), C(P), TAL(CATTOT)

       INTEGER H, HCUM, J
       REAL DIF

************************************************************************
*           Begin.
************************************************************************
       DSM = 0.
       HCUM = 0
       DO J = 1, P
         IF (CAT(J) .EQ. 0) THEN
           IF (T(J) < LIT) THEN                ! scalar
             DIF = X(J) - C(J)
           ELSE                                ! cyclic
             DIF = MOD(ABS(X(J) - C(J)), T(J))
             DIF = MIN(DIF, T(J) - DIF)              ! choose the shorter path
           END IF
         ELSE
           H = HCUM + NINT(X(J))
           IF (T(J) > -2.5) THEN                  ! categorical
             DIF = (POP - TAL(H)) / POP
           ELSE                             ! contra-categorical
             DIF = (TAL(H) - SELF) / POP
           END IF
           HCUM = HCUM + CAT(J)
         END IF  ! cat > 0 ?
         IF (ROBUST) THEN
           DSM = DSM + ABS(W(J) * DIF)     ! Manhattan dissimilarity
         ELSE
           DSM = DSM + (W(J) * DIF) **2    ! Squared Euclidean dissimilarity
         END IF
       
         IF (DSM > 1E30) THEN 
           PRINT *, 'dsm: huge value.  X: ', X, ' p: ', P, ' T: ', T,
     $        ' W: ', W, ' cat: ', CAT, ' Cattot: ', CATTOT, ' tal:',
     $       TAL, ' pop:', POP, ' lit: ', LIT, ' self:', SELF,
     $       ' robust: ', ROBUST, ' best: ', BEST
         END IF  
       
         IF (DSM .GE. BEST) RETURN     ! bail out
       END DO  ! next j
       RETURN
      END  ! of DSM


************************************************************************
* theur - threshhold heuristic, the maximin dissimilarity
*
*___Name____________Type______I/O_______Description________________________
*   X(PLIM,N)       Real      In        Data points
*   PLIM            Integer   In        Row dimension of X
*   P               Integer   In        Dimension of the data
*   N               Integer   In        Number of points
*   T(P)            Real      In        Periods of cyclic data
*   W(P)            Real      In        Importance weights
*   CAT(P)          Integer   In        # categories each variable
*   CATTOT          Integer   In        Total categories
*   C(PLIM,K)       Real      Neither   Center points of clusters
*   K               Integer   In        Default number of clusters
*   TAL(CATTOT,K)   Real      Neither   Observer frequencies
*   Z(N)            Integer   Neither   Argument to CENTER
*   POP(K)          Real      Neither   Population each cluster
*   ROBUST          Logical   In        Use medians instead of means
*   LIVEC(K)        Integer   Neither   Argument to CENTER
*   IPERM(N)        Integer   Neither   Permutation vector
*   SHORTX(N)       Real      Neither   Dissimilarity to each point
*   WORK(2*N)       Real      Neither   Argument to CENTER
*   THRESH          Real      Out       Threshhold distance
************************************************************************
      SUBROUTINE THEUR (X, PLIM, P, N, T, W, CAT, CATTOT, C, K, TAL, Z, 
     $  POP, ROBUST, LIVEC, IPERM, SHORTX, WORK, THRESH)
       IMPLICIT NONE
       INTEGER PLIM, P, N, CAT, CATTOT, K, Z, LIVEC, IPERM
       REAL X, T, W, C, TAL, POP, SHORTX, WORK, THRESH
       LOGICAL ROBUST
       DIMENSION X(PLIM,N), T(P), W(P), CAT(P), C(PLIM,K), 
     $  TAL(CATTOT,K), Z(N), POP(K), LIVEC(K), IPERM(N), SHORTX(N), 
     $  WORK(2*N)

*          Local variables
       INTEGER H, HCUM, I, I0, I1, J, L
       REAL BIG, D, FAR, LIT, SELF

*          External functions
       REAL DSM, RMACHINE

       DATA I1 / 1 /  ! avoid compiler warning
       
************************************************************************
*            Begin.
************************************************************************
       LIT = RMACHINE(2)
       BIG = RMACHINE(3)
       DO I = 1, N
         SHORTX(I) = BIG
       END DO

*         First center is the overall average
       DO I = 1, N
         Z(I) = 1
       END DO
       LIVEC(1) = 1
       CALL CENTER (X, PLIM, P, N, T, CAT, CATTOT, C, 1, TAL, Z, POP, 
     $  ROBUST, .FALSE., LIVEC, IPERM, WORK)
       I0 = 0    ! remember what point was chosen

*        Choose each successive center as the point with the longest
*         minimum distance to an existing center
       DO L = 2, K
         FAR = 0.
         I0 = I1
         DO I = 1, N            ! compare each point to previous center
           IF (I .EQ. I0) THEN
             SELF = 1.
           ELSE
             SELF = 0.
           END IF
           D = DSM (X(1,I), P, T, W, CAT, CATTOT, C(1,L-1),
     $             TAL(1,L-1), POP(L-1), LIT, SELF, ROBUST, SHORTX(I))
           IF (D < SHORTX(I)) SHORTX(I) = D
           IF (SHORTX(I) > FAR) THEN           ! longest short distance
             FAR = SHORTX(I)
             I1 = I
           END IF
         END DO                      ! next data point
         
         DO H = 1, CATTOT
           TAL(H,L) = 0.
         END DO
         HCUM = 0.
         DO J = 1, P         ! assign center
           C(J,L) = X(J,I1)
           IF (CAT(J) > 0) THEN
             H = HCUM + NINT(X(J,I1))
             TAL(H,L) = 1.
             HCUM = HCUM + CAT(J)
           END IF
         END DO
         POP(L) = 1.
       END DO               ! next center to initialize

*   Threshhold distance is longest distance of any point to its center
       THRESH = 0.
       DO I = 1, N
         IF (I .EQ. I1) THEN
           SELF = 1.
         ELSE
           SELF = 0.
         END IF
         D = DSM (X(1,I), P, T, W, CAT, CATTOT, C(1,K), TAL(1,K),
     $     POP(K), LIT, SELF, ROBUST, SHORTX(I))
         IF (D < SHORTX(I)) SHORTX(I) = D
         IF (SHORTX(I) > THRESH) THRESH = SHORTX(I)
       END DO
       RETURN
      END  ! of THEUR


************************************************************************
* Initialize cluster centers based on distance
*
* Refer to:
*    "kmeans++: the advantages of careful seeding"
*    David Arthur and Sergei Vassilvitskii
*    Proceedings of the eighteenth annual ACM-SIAM Symposium
*      on Discrete Algorithms, 2007
*
*___Name____________Type______I/O_______Description_____________________
*   X(PLIM,N)       Real      In        Data points
*   PLIM            Integer   In        Row dimension of X
*   P               Integer   In        Dimension of the data
*   N               Integer   In        Number of points
*   T(P)            Real      In        Period of modulo arithmetic
*   W(P)            Real      In        Importance weights
*   CAT(P)          Integer   In        Categories each variable
*   CATTOT          Integer   In        Total categories
*   C(PLIM,K)       Real      Out       Center points of clusters
*   K               Integer   In        Number of clusters to make
*   TAL(CATTOT,K)   Real      Out       Category frequencies each cluster
*   POP(K)          Real      Out       Populations 
*   ROBUST          Logical   In        Use medians instead of means
*   SHORTX(N)       Real      Neither   Squared distance to each point
*   CUMD(N)         Real      Neither   Cumulative squared distance
*   THRESH          Real      Out       Threshhold distance
************************************************************************
      SUBROUTINE DINIT (X, PLIM, P, N, T, W, CAT, CATTOT, C, K, TAL, 
     $  POP, ROBUST, SHORTX, CUMD, THRESH)

       IMPLICIT NONE
       INTEGER PLIM, P, N, CAT, CATTOT, K
       REAL X, T, W, C, TAL, POP, SHORTX, CUMD, THRESH
       LOGICAL ROBUST
       DIMENSION X(PLIM,N), T(P), W(P), CAT(P), C(PLIM,K), 
     $  TAL(CATTOT,K), POP(K), SHORTX(N), CUMD(N)

*          Local variables
       INTEGER H, HCUM, I, I1, J, L
       REAL BIG, D, LIT, SELF, TOT, U

*          External functions
       INTEGER HAT
       REAL DSM, 
     $      RMACHINE,
     $      URAND

       LIT = RMACHINE(2)
       BIG = RMACHINE(3)

* Begin.  Choose each center with probability proportional to its
* least squared distance from any existing center.
       DO I = 1, N
         SHORTX(I) = BIG
       END DO

       I1 = HAT(N)  ! choose first center at random
       DO H = 1, CATTOT
         TAL(H,1) = 0.
       END DO
       HCUM = 0
       DO J = 1, P                         ! copy
         C(J,1) = X(J,I1)
         IF (CAT(J) > 0) THEN
           H = HCUM + NINT(X(J,I1))
           TAL(H,1) = 1.
           HCUM = HCUM + CAT(J)
         END IF
       END DO
       POP(1) = 1.

       DO L = 2, K
         TOT = 0.
         DO I = 1, N             ! compare each point to previous center
           IF (I .EQ. I1) THEN
             SELF = 1.
           ELSE
             SELF = 0.
           END IF
           D = DSM (X(1,I), P, T, W, CAT, CATTOT, C(1,L-1), TAL(1,L-1),
     $                POP(L-1), LIT, SELF, ROBUST, SHORTX(I))
     
            IF (D > 1E30) PRINT *, 'dinit: huge D!!!'
     
           IF (D < SHORTX(I)) SHORTX(I) = D
           TOT = TOT + SHORTX(I)
           CUMD(I) = TOT
         END DO                      ! next data point

         U = URAND() * TOT    ! uniform at random over cumulative distance
         CALL RBSGT (CUMD, N, U, I1)
         DO H = 1, CATTOT
           TAL(H,L) = 0.
         END DO
         HCUM = 0
         DO J = 1, P         ! assign center
           C(J,L) = X(J,I1)
           IF (CAT(J) > 0) THEN
             H = HCUM + NINT(X(J,I1))
             TAL(H,L) = 1.
             HCUM = HCUM + CAT(J)
           END IF
         END DO
         POP(L) = 1.
       END DO               ! next center to initialize

*   Threshhold distance is longest distance of any point to its center
       THRESH = 0.
       DO I = 1, N
         IF (I .EQ. I1) THEN
           SELF = 1.
         ELSE
           SELF = 0.
         END IF    
         D = DSM (X(1,I), P, T, W, CAT, CATTOT, C(1,K), TAL(1,K),
     $                POP(K), LIT, SELF, ROBUST, SHORTX(I))
         IF (D < SHORTX(I)) SHORTX(I) = D
         IF (SHORTX(I) > THRESH) THRESH = SHORTX(I)
       END DO
       
       RETURN
      END  ! of DINIT



***********************************************************************
* KCLUST - partition clustering for mixed data
* origin:  Hugo Steinhaus, 1956
*
*___Variable___________Type______I/O_______Description__________________
*   X(PLIM,N)          Real      In        Data points
*   PLIM               Integer   In        Row dimension of x
*   P                  Integer   In        Dimension of the data
*   N                  Integer   In        Number of points
*   T(P)               Real      In        Period of modulo arithmetic
*   W(P)               Real      In        Importance weights
*   CAT(P)             Integer   In        # categories each variable
*   CATTOT             Integer   In        Total categories
*   C(PLIM,KMAX)       Real      Both      Center points of clusters
*   KMIN               Integer   In        Fewest clusters to make
*   K                  Integer   Both      Number of clusters to make
*   KMAX               Integer   In        Most clusters to make
*   TAL(CATTOT,KMAX)   Real      Both      Observer frequencies
*   Z(N)               Integer   Both      What cluster a point is in
*   POP(KMAX)          Real      Both      # of points in each cluster
*   ROBUST             Logical   In        Use medians
*   THRESH             Real      In        Threshhold for forming new cluster
*   LIVEC(KMAX)        Integer   Neither   Workspace, 1 = live, 0 = dead
*   LIVEX(N)           Integer   Neither   Workspace, time-to-live each point
*   IPERM(N)           Integer   Neither   Permutation vector
*   SHORTX(N)          Real      Neither   Shortest distance to current center
*   WORK(2*N)          Real      Neither   Workspace for averages
*   U                  Real      Out       Residual
*   IFAULT             Integer   Out       Error code
*                                            0 = no errors
*                                            1 = too few clusters to permit 
*                                                 a cluster to be deleted
*                                            2 = too many iterations
*                                            3 = K out of bounds
*                                            4 = can't measure point
*                                            5 = can't insert new cluster
************************************************************************
      SUBROUTINE KCLUST (X,PLIM,P,N,T,W,CAT,CATTOT,C,KMIN,K,KMAX,
     & TAL,Z,POP,ROBUST,THRESH,LIVEC,LIVEX,IPERM,SHORTX,WORK,U,IFAULT)

       IMPLICIT NONE
       INTEGER PLIM, P, N, CAT, CATTOT, KMIN, K, KMAX, Z,
     &          LIVEC, LIVEX, IPERM, IFAULT
       REAL X, T, W, C, TAL, POP, THRESH, SHORTX, WORK, U
       LOGICAL ROBUST
       DIMENSION X(PLIM,N), T(P), W(P), CAT(P), C(PLIM,KMAX),
     &   TAL(CATTOT,KMAX), Z(N), POP(KMAX),
     &   LIVEC(KMAX), LIVEX(N), IPERM(N), SHORTX(N), WORK(2*N)

*               Constants
       INTEGER HS, ITER
       PARAMETER (HS = 13,        ! length of history hash buffer
     &            ITER =  1000)   ! maximum iterations

*                Local variables
       INTEGER
     &         G, H, HCUM, I, I1, IINS, J, L,   ! counters
     &         HL,                              ! length in bytes of Z[]
     &         K0,                              ! temp. value of K
     &         LOLD, LNEW                       ! present, new cluster ID

       INTEGER*4 HASH, HIST(HS)                 ! hashes of prior Z()

       REAL BEST, WORST,    ! Shortest, longest distance to a center
     &      BIG, LIT,       ! very large number, small #
     &      D,              ! a distance
     &      SELF            ! this cluster already contains this point

       LOGICAL CHANGE             ! whether any points have been reassigned

*              external functions
        INTEGER BITSZ                       ! replacement for BIT_SIZE
        REAL DSM, RMACHINE
        LOGICAL ALMEQ

        DATA IINS / 1 /    ! avoid compiler warning

************************************************************************
*           Begin
************************************************************************
       IFAULT = 0                          ! validate, initialize
       LIT = RMACHINE(2)
       BIG = RMACHINE(3)
       IF (K < 1 .OR. K > N) THEN
         IFAULT = 3
         PRINT *, 'kclust:  k= ', K, ' but n=', N
         RETURN
       END IF
       HL = N * BITSZ() / 8

*          initialize
       IF (K .EQ. 1) GO TO 90
       DO L = 1, K                  ! initialize live sets
         LIVEC(L) = 2
       END DO
       DO I = 1, N
         LIVEX(I) = 1
         SHORTX(I) = BIG
       END DO

************************************************************************
*                      Main loop
************************************************************************
       DO G = 1, ITER
         CHANGE = .FALSE.
         CALL SHUFFL (IPERM, N)

************************************************************************
*            Find nearest center for each point
************************************************************************
         WORST = 0.
         DO I1 = 1, N
           I = IPERM(I1)
           LOLD = Z(I)                ! old assignment
           LNEW = LOLD
           BEST = SHORTX(I)

           DO 30 L = 1, K
             IF (LIVEX(I) .LE. 0 .AND. LIVEC(L) .LE. 0) GO TO 30
             IF (L .EQ. LOLD) THEN
               SELF = 1.
             ELSE
               SELF = 0.
             END IF
             D = DSM (X(1,I), P, T, W, CAT, CATTOT, C(1,L),
     $                TAL(1,L), POP(L), LIT, SELF, ROBUST, BEST)

             IF (D < BEST) THEN          ! new nearest center
               BEST = D
               LNEW = L
             END IF
  30       CONTINUE        ! next L
           IF (ALMEQ(BEST, BIG)) THEN   ! incomparable point?
             IFAULT = 4
             RETURN
           END IF
           SHORTX(I) = BEST
           
           IF (BEST > WORST) THEN  ! Compare distance to new cluster threshhold
             IINS = I
             WORST = BEST
           END IF           

************************************************************************
*              Reassign point
************************************************************************
           IF (LOLD .NE. LNEW) THEN
             Z(I) = LNEW
             LIVEC(LOLD) = 2           ! recenter old cluster
             LIVEC(LNEW) = 2           ! recenter new cluster
             CHANGE = .TRUE.
           END IF
         END DO         ! next i1

************************************************************************
* If shortest distance exceeds threshhold, make a new cluster
************************************************************************
         IF (WORST > THRESH .AND. K < KMAX) THEN
           IF (K .GE. KMAX) THEN
             IFAULT = 5
             RETURN
           END IF
           K = K + 1
           LOLD = Z(IINS)
           LNEW = K
           DO J = 1, P
             C(J,LNEW) = X(J,IINS)
           END DO
           DO H = 1, CATTOT
             TAL(H,LNEW) = 0.
           END DO
           HCUM = 0
           DO J = 1, P
             IF (CAT(J) > 0) THEN
               H = HCUM + NINT(X(J,IINS))
               TAL(H,LNEW) = 1.
               HCUM = HCUM + CAT(J)
             END IF
           END DO
           Z(IINS) = LNEW
           LIVEC(LOLD) = 2
           LIVEC(LNEW) = 2
           CHANGE = .TRUE.
         END IF

         IF (.NOT. CHANGE) GO TO 90       ! algorithm halts

************************************************************************
*           Find cluster centers
************************************************************************
         DO L = 1, K                            ! decrement LIVEC
           LIVEC(L) = LIVEC(L) - 1
         END DO
         CALL CENTER (X, PLIM, P, N, T, CAT, CATTOT,
     &    C, K, TAL, Z, POP, ROBUST, .FALSE., LIVEC, IPERM, WORK)

************************************************************************
*          Delete empty clusters
************************************************************************
         K0 = K
         DO L = K0, 1, -1
           IF (POP(L) < 0.5) THEN
             IF (K .LE. KMIN) THEN
               IFAULT = 1
               RETURN
             END IF
             DO J = 1, P
               C(J,L) = C(J,K)
             END DO
             DO H = 1, CATTOT
               TAL(H,L) = TAL(H,K)
             END DO
             DO I = 1, N
               IF (Z(I) .EQ. K) Z(I) = L
              END DO
             POP(L) = POP(K)
             LIVEC(L) = LIVEC(K)
             POP(K) = 0.
             K = K - 1
           END IF ! empty cluster?
         END DO  ! next L

************************************************************************
*            Check if anything is changing
************************************************************************
         HASH = -2128831035           ! canonical arbitrary initialization
         CALL FNV32(Z, HL, HASH)      ! hash of membership vector
         DO J = 1, MIN(G-1, HS)
           IF (HASH .EQ. HIST(J)) GO TO 90      ! exit on loop criterion
         END DO
         J = MOD(G, HS)
         IF (J .EQ. 0) J = HS
         HIST(J) = HASH                          ! remember

************************************************************************
*          Determine live points
************************************************************************
        DO I = 1, N
          LIVEX(I) = 0
        END DO
        SELF = 1.
        DO I = 1, N
          L = Z(I)
          IF (LIVEC(L) > 0) THEN
            D = DSM (X(1,I), P, T, W, CAT, CATTOT, C(1,L),
     $               TAL(1,L), POP(L), LIT, SELF, ROBUST, BIG)
            IF (D > SHORTX(I)) LIVEX(I) = 1
            SHORTX(I) = D
          END IF
        END DO

       END DO                   ! next g
       IFAULT = 2               ! too many iterations

************************************************************************
*         Compute residual
************************************************************************
  90   U = 0.
       SELF = 1.
       DO I = 1, N
         L = Z(I)
         D = DSM (X(1,I), P, T, W, CAT, CATTOT, C(1,L),
     $               TAL(1,L), POP(L), LIT, SELF, ROBUST, BIG)
         SHORTX(I) = D
         U = U + D
       END DO  ! next i
       RETURN
      END  ! of kclust



*--------------------- --------------------------------------------------
* jclust - precedes kclust
*
*___Variable___________Type______I/O______Description___________________
*   X(PLIM,N)          Real      In       Data points in R^p
*   PLIM               Integer   In       Row dimension of X
*   P                  Integer   In       Measurements per point
*   N                  Integer   In       Number of points
*   T(P)               Real      In       Period of modulo arithmetic
*   W(P)               Real      In       Importance weights of variables
*   CAT(P)             Integer   In       Number of categories
*   CATTOT             Integer   In       The total of CAT(1:p)
*   M                  Integer   In       Number of sources
*   C(PLIM,KMAX)       Real      Out      Cluster centroids
*   KMIN               Integer   In       Fewest clusters to make
*   K                  Integer   In/Out   Default/actual clusters to make
*   KMAX               Integer   In       Most clusters to make
*   TAL(CATTOT,KMAX)   Real      Work     Observer frequencies
*   Z(N)               Integer   Out      What cluster each point is in
*   POP(KMAX)          Real      Out      Number of points in each cluster
*   OPT                Integer   In       Options
*   ALG                Integer   In       Select algorithm {0,1}
*   ZSAVE(N,CTRIES)    Integer   Work     Best Z
*   IPERM(N)           Integer   Work     Permutation vector
*   LIVEX(N)           Integer   Work     Live set membership, each point
*   LIVEC(KMAX)        Integer   Work     Live set membership, each cluster
*   SHORTX(N)          Real      Work     Distance to center, each point
*   XW(2*N)            Real      Work     Work array for averages
*   IERR               Integer   Out      Error code
*                                           0 = no errors
*                                           {1,2,3,4,5} = see KCLUST
*                                           6 = K out of bounds
*                                           7 N/A (see package CLUELA)
*                                           8 = bad array bounds
*                                           9 = no weights assigned
*                                          10 N/A (see package CLUELA)
*                                          11 = negative weights
*                                          12 = bad value in CAT[]
*                                          13 = tclust fails
*                                          67 = bad value for ALG
*-----------------------------------------------------------------------
      SUBROUTINE JCLUST (X, PLIM, P, N, T, W, CAT, CATTOT, M, 
     $  C, KMIN, K, KMAX, TAL, Z, POP, OPT, ALG, ZSAVE, IPERM,
     $  LIVEX, LIVEC, SHORTX, XW, IERR)

       IMPLICIT NONE
       INTEGER CTRIES
       PARAMETER (CTRIES = 10)    ! # repeats of k-means
       INTEGER PLIM, P, N, CAT, CATTOT, M, KMIN, K, KMAX, Z,
     $  OPT, ALG, ZSAVE, IPERM, LIVEX, LIVEC, IERR
       REAL X, T, W, C, TAL, POP, SHORTX, XW
       DIMENSION X(PLIM,N), T(P), W(P), CAT(P), C(PLIM,KMAX),
     $  TAL(CATTOT,KMAX), Z(N), POP(KMAX), ZSAVE(N,CTRIES), IPERM(N),
     $  LIVEX(N), LIVEC(KMAX), SHORTX(N), XW(2*N)

*-----------------------------------------------------------------------
*          Compiled-in options
*-----------------------------------------------------------------------
       INTEGER TTRIES
       LOGICAL CRISP, REEST
       PARAMETER (CRISP = .FALSE.,       ! use all-or-nothing centers
     $            REEST = .FALSE.,        ! change input array
     $            TTRIES = 5)

*-----------------------------------------------------------------------
*            local variables
*-----------------------------------------------------------------------
      INTEGER
     $        H, i,             ! count trials, count points
     $        HBEST,            ! best trial
     $        HOK,              ! successful trials
     $        J,                ! count dimensions
     $        JTIES,            ! ties in IMODE
     $        KBEST, KTRY,      ! possible value for k
     $        KARR(CTRIES),     ! results for k, each trial
     $        KAVE(CTRIES),     ! workspace for IMODE
     $        KPERM(CTRIES),    ! workspace for IMODE
     $        L                 ! count clusters

      REAL
     $     BIG,       ! very large number
     $     FM,              ! M as Real
     $     LIT,             ! very small number
     $     D,               ! a dissimilarity
     $     SELF,            ! argument for DSM
     $     THRESH, T0,      ! threshhold for splitting
     $     U,               ! residual from clustering routines
     $     UARR(CTRIES)     

       LOGICAL
     $        BRIEF,      ! skip clustering restarts
     $        HYP,        ! hyperplane initialization
     $        ROBUST      ! robust statistics

*             external functions
       INTEGER HAT              ! random integer 1...n
       REAL    DIST1, DIST2, 
     $         DSM,   ! distance functions
     $         RMACHINE             ! find machine constants
!!       LOGICAL SAFDIV
       EXTERNAL DIST1, DIST2

*-----------------------------------------------------------------------
*       begin executable statements
*-----------------------------------------------------------------------
       IERR = 0                       ! clear error codes
       LIT = RMACHINE(2)              ! very small number
       BIG = RMACHINE(3)              ! very large number
       FM = FLOAT(M)
       KBEST = K                      ! avoid compiler warning

*-----------------------------------------------------------------------
*            validate
*-----------------------------------------------------------------------
       IF (KMIN < 1 .OR. K < KMIN .OR. K > KMAX .OR. KMAX > N) THEN
         PRINT *, 'jclust!! kmin, k, kmax ', KMIN, K, KMAX, ' n=', N
         IERR = 6
         RETURN
       END IF

       IF (P < 1 .OR. N < 1) THEN
         PRINT *, 'jclust! p: ', P, ' n: ', N
         IERR = 8
         RETURN
       END IF
              
*        at least one element of W[] must be nonzero
       DO J = 1, P
         IF (W(J) > LIT) GO TO 10
       END DO
       IERR = 9
       RETURN

  10   DO J = 1, P
         IF (W(J) < 0.) THEN
           PRINT *, 'jclust!! negative weights! W: ', W
           IERR = 11
           RETURN
         END IF
       END DO

*    Ensure that categorical and contra-categorical data are positive integers
       DO i = 1, N
         DO J = 1, P
           IF (T(J) < -0.5) THEN
             IF (X(J,i) < 0.5 .OR. X(J,i) > FLOAT(CAT(J)) + 0.4999) THEN
               PRINT *, 'CAT(', J, ')=', CAT(J), 
     $          ' but X(', J, ',', I, ')=', X(J,i)
               IERR = 12
               RETURN
             END IF
           END IF
         END DO
       END DO       

*-----------------------------------------------------------------------
*           decode input
*-----------------------------------------------------------------------
       ROBUST = BTEST(OPT, 0)       ! medians instead of means
       BRIEF = BTEST(OPT, 1)        ! skip clustering restarts
       HYP = BTEST(OPT, 2)

*-----------------------------------------------------------------------
*        special case k = n
*-----------------------------------------------------------------------
       IF (KMIN .EQ. N) THEN
         DO i = 1, N
           DO J = 1, P
             C(J,i) = X(J,i)
           END DO
         END DO
         DO i = 1, N
           Z(i) = i
         END DO
         DO L = 1, K
           POP(L) = 1.
         END DO
         RETURN
       END IF

*-----------------------------------------------------------------------
*         special case k = 1
*-----------------------------------------------------------------------
       IF (KMAX .EQ. 1) THEN
         DO i = 1, N                  ! all memberships the same
           Z(i) = 1
         END DO
         POP(1) = N
         GO TO 90
       END IF

*-----------------------------------------------------------------------
*        choose what algorithm
*-----------------------------------------------------------------------
       IF (ALG .EQ. 0 .OR. ALG .EQ. 1) THEN
         CONTINUE
       ELSE
         PRINT *, 'jclust: bad value for alg!'
         IERR = 67
         RETURN
       END IF

*-----------------------------------------------------------------------
*         initialize k-means routines, CTRIES attempts, keep the best
*-----------------------------------------------------------------------
       HOK = 0
       DO H = 1, CTRIES            ! seek least residual
         IF (HYP) THEN
           IF (ROBUST) THEN
             CALL HCLUST (X,PLIM,P,N,T,W,K,Z,POP,IPERM,LIVEC,DIST1,IERR)
           ELSE
             CALL HCLUST (X,PLIM,P,N,T,W,K,Z,POP,IPERM,LIVEC,DIST2,IERR)
           END IF
           IF (IERR > 0) PRINT *, 'hclust: Trouble partitioning'

*            find initial cluster centers
           DO L = 1, K               ! all clusters live
             LIVEC(L) = 1
           END DO
           CALL CENTER (X, PLIM, P, N, T, CAT, CATTOT,
     $         C, K, TAL, Z, POP, ROBUST, .FALSE., LIVEC, IPERM, XW)
           THRESH = 0.
           SELF = 1.
           DO I = 1, N
             L = Z(I)
             D = DSM (X(1,I), P, T, W, CAT, CATTOT, C(1,L), TAL(1,L),
     $               POP(L), LIT, SELF, ROBUST, BIG)
             IF (D > THRESH) THRESH = D
           END DO
         ELSE
           IF (H .EQ. 1 .OR. REEST) THEN
             CALL THEUR (X, PLIM, P, N, T, W, CAT, CATTOT, C, K, TAL, Z,
     $             POP, ROBUST, LIVEC, IPERM, SHORTX, XW, THRESH)
!            PRINT *, 'theur returns!  thresh= ', THRESH
           END IF
     
           CALL DINIT (X, PLIM, P, N, T, W, CAT, CATTOT, C, K, TAL,
     $       POP, ROBUST, SHORTX, XW, T0)
!           PRINT *, 'dinit returns.  t0 = ', T0
           DO i = 1, N
             Z(i) = HAT(K)
           END DO
         END IF

!!!!!!!!!!!!!!!
!         PRINT *, 'calculated thresh as: ', THRESH
!                use thresh determined by DINIT ?
!!         IF (.NOT. HYP .AND. ALG .EQ. 0) THRESH = T0

*-----------------------------------------------------------------------
*           Cluster
*-----------------------------------------------------------------------
         KTRY = K
         U = BIG
         CALL KCLUST (X,PLIM,P,N,T,W,CAT,CATTOT,C,KMIN,KTRY,KMAX,
     $      TAL,Z,POP,ROBUST,THRESH,LIVEC,LIVEX,IPERM,SHORTX,XW,U,IERR)

!!!!!       debug:
!           PRINT *, 'returned from kclust.  h=', h, ' u=', u, ' k=',
!     $         KTRY, ' ierr=', IERR

         IF (BRIEF .AND. IERR .EQ. 0) THEN
           KBEST = KTRY
           GO TO 50      ! break from loop H
         END IF         

         IF (IERR .EQ. 0) THEN
           HOK = HOK + 1
           KARR(HOK) = KTRY
           UARR(HOK) = U
           DO i = 1, N
             ZSAVE(i,HOK) = Z(i)
           END DO
           
*-----------------------------------------------------------------------
* Adjust X towards an approximation of C as a function of X. This step is 
* intended to have the same effect as Dynamic Time Warping, since 
* data clustering by itself does not take advantage of information about 
* proximity of events in their original context.
*----------------------------------------------------------------------
           IF (REEST .AND. ALG .EQ. 0) THEN
             DO i = 1, N                      ! weight by inverse distance
               SHORTX(i) = FM / (FM + SHORTX(i))
             END DO
             CALL RECS (X, PLIM, P, N, T, M, C, KTRY, Z, IPERM, LIVEX,
     $                   SHORTX, XW(1), XW(N+1))
           END IF  ! re-estimate?

         ELSE
           CONTINUE
           PRINT *, 'kclust:  error#', IERR
         END IF  ! successful kclust run?
       END DO  ! next H

       IF (HOK .EQ. 0) RETURN  ! failure
       CALL IMODE (KARR, KAVE, KPERM, HOK, JTIES)  ! most common K
       U = BIG
       HBEST = 1

*         choose lowest residual from trials with modal KTRY
       IF (JTIES .EQ. 1) THEN
         KBEST = KAVE(1)
         DO H = 1, HOK
           IF (KARR(H) .EQ. KBEST .AND. UARR(H) < U) THEN
             HBEST = H
             U = UARR(H)
           END IF
         END DO

*         choose lowest residual from any trial
       ELSE
         DO H = 1, HOK
           IF (UARR(H) < U) THEN
             HBEST = H
             KBEST = KARR(H)
             U = UARR(H)
           END IF
         END DO
       END IF

*          restore best saved values
       PRINT *, 'choosing results of trial ', HBEST
       DO i = 1, N
         Z(i) = ZSAVE(i,HBEST)         
       END DO

*-----------------------------------------------------------------------
*           Come here when skipping restarts
*-----------------------------------------------------------------------
  50   IERR = 0
       K = KBEST
       IF (ALG .EQ. 0 .OR. K .EQ. 1 .OR. K .EQ. N) GO TO 90

*-----------------------------------------------------------------------
*              Mixture decomposition.
*-----------------------------------------------------------------------
       DO H = 1, TTRIES
         KTRY = K
         IF (REEST .AND. H < TTRIES) THEN
           J = 3
         ELSE
           J = 500
         END IF
         CALL TCLUST (X, PLIM, P, N, T, CAT, CATTOT, M, C, KMIN, KTRY,
     $  KMAX, TAL, Z, POP, ROBUST, IPERM, XW, XW(N+1), SHORTX, J, IERR)
         PRINT *, 'jclust iter ', H, ' tclust returns #', IERR
         IF (IERR .NE. 0) THEN
!!!         PRINT *, 'jclust: returning failure!!'
           IERR = 13
!!         PAUSE 'tclust has failed!'
!          RETURN     ! failure
         END IF

         IF (IERR .EQ. 0 .AND. REEST .AND. H < TTRIES) THEN
!           DO L = 1, KTRY           ! Crispify
!             LIVEC(L) = 1
!           END DO
!           CALL CENTER (X, PLIM, P, N, T, CAT, CATTOT, 
!     $          C, KTRY, TAL, Z, POP, ROBUST, .FALSE., LIVEC, IPERM, XW)

           SELF = 1.
           DO i = 1, N                      ! weight by inverse distance
             L = Z(i)
             D = DSM (X(1,I), P, T, W, CAT, CATTOT, C(1,L), TAL(1,L),
     $             POP(L), LIT, SELF, ROBUST, BIG)
             SHORTX(i) = FM / (FM + D)
           END DO
           CALL RECS (X, PLIM, P, N, T, M, C, KTRY, Z, IPERM, LIVEX,
     $                   SHORTX, XW(1), XW(N+1))
         ELSE IF (IERR .NE. 0) THEN
           CONTINUE
         ELSE
           GO TO 70
         END IF  ! re-estimate?
       END DO

  70   K = KTRY
       IF (CRISP) THEN
         CONTINUE
       ELSE
         RETURN   ! success
       END IF

*-----------------------------------------------------------------------
*         Count cluster populations, return centers
*-----------------------------------------------------------------------
  90   DO L = 1, K         ! all clusters live
         LIVEC(L) = 1
       END DO
       CALL CENTER (X, PLIM, P, N, T, CAT, CATTOT,
     $    C, K, TAL, Z, POP, ROBUST, .TRUE., LIVEC, IPERM, XW)
       RETURN
      END  ! of JCLUST
*++++++++++++++++++++++++ End of file cluster.f ++++++++++++++++++++++++
