*------------------------------------------------------------------------
* The Common Algorithm for data clustering
* by Andy Allinger, 2011-2017, released to the public domain
* This program may be used by any person for any purpose.
*
* Refer to:
*   Robust Cluster Analysis Via Mixtures Of Multivariate t-Distributions
*   Geoffrey J. McLachlan and Donald Peel
*   University of Queensland
*
* and:
*   Roubust mixture modelling using the t distribution
*   D. Peel and G. J. McLachlan
*   Statistics and Computing v.10 pp.339-348, 2000.
*------------------------------------------------------------------------

*-----------------------------------------------------------------------
* dtpfc - distance to point from center
*
*___Name__________Type_______________In/Out___Description_______________
*   X(P)          Real               In       Datum vector
*   P             Integer            In       Length of X
*   T(P)          Real               In       Period of modulo arithmetic
*   CAT(P)        Integer            In       # categories each variable
*   CATTOT        Integer            In       Total number of categories
*   LIMCON        Integer            In       # continuous variables
*   LIMCAT        Integer            In       # categorical variables
*   C(P)          Real               In       Centroid vector
*   TAL(CATTOT)   Real               In       Count of categorical data
*   POP           Real               In       Cluster population
*   LIT           Real               In       Negligible value
*   DS(LIMCON)    Double Precision   Out      Difference, scalars
*   DT(LIMCAT)    Double Precision   Out      Difference, categoricals
*-----------------------------------------------------------------------
      SUBROUTINE DTPFC (X, P, T, CAT, CATTOT, LIMCON, LIMCAT, 
     $                    C, TAL, POP, LIT, DS, DT)
       IMPLICIT NONE
       INTEGER P, CAT, CATTOT, LIMCON, LIMCAT
       REAL X, T, C, TAL, POP, LIT
       DOUBLE PRECISION DS, DT
       DIMENSION X(P), T(P), CAT(P), C(P), TAL(CATTOT),
     $             DS(LIMCON), DT(LIMCAT)

*            local variables
       INTEGER H, HCUM, J, JS, JT
       DOUBLE PRECISION DIF, HT
       
*-----------------------------------------------------------------------
*          begin
*-----------------------------------------------------------------------
       HCUM = 0
       JS = 0
       JT = 0
       DO J = 1, P
         IF (CAT(J) .EQ. 0) THEN
           JS = JS + 1
           IF (T(J) < LIT) THEN                ! scalar
             DS(JS) = X(J) - C(J)
           ELSE                                ! modulo
             HT = 0.5 * T(J)
             DIF = MOD(X(J) - C(J), T(J))
             IF (DIF < -HT) THEN
               DIF = DIF + T(J)
             ELSE IF (DIF > HT) THEN
               DIF = DIF - T(J)
             END IF
             DS(JS) = DIF
           END IF
         ELSE
           JT = JT + 1
           H = HCUM + NINT(X(J))
           IF (T(J) > -2.5) THEN             ! categorical
             DIF = (TAL(H) + 1.) / (POP + FLOAT(CAT(J)))  ! Laplace smoothing
           ELSE                              ! contra-categorical
             DIF = (POP - TAL(H) + 1.) / 
     $               ((POP + 1.) * FLOAT(CAT(J)) - POP)
           END IF
           DT(JT) = DIF
           HCUM = HCUM + CAT(J)
         END IF ! cat > 0 ?
       END DO ! next j
       RETURN
      END ! of DTPFC


************************************************************************
* tclust - Mixture decomposition with the Student t distribution
*
*___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
*   CAT(P)             Integer   in     # categories each variable 
*   CATTOT             Integer   in     total # categories
*   M                  Integer   in     number of input sources
*   C(PLIM,KMAX)       Real      out    cluster centers
*   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      work   category counts each cluster
*   Z(N)               Integer   out    maximal cluster membership each point
*   POP(KMAX)          Real      work   cluster population
*   ROBUST             Logical   in     use weighted median
*   IND(N)             Integer   in     permutation vector
*   WORK(N)            Real      work   temp space of X
*   QW(N)              Real      work   weights of weighted averages
*   AVE(N)             Real      work   result array for WCYMEN
*   ITER               Integer   in     quit after this many iterations
*   IERR               Integer   out    error code
************************************************************************
      SUBROUTINE TCLUST (X, PLIM, P, N, T, CAT, CATTOT, M, C, KMIN, K, 
     $  KMAX, TAL, Z, POP, ROBUST, IND, WORK, QW, AVE, ITER, IERR)

       IMPLICIT NONE
       INTEGER PLIM, P, N, CAT, CATTOT, M, KMIN, K, KMAX, Z, IND, ITER, 
     $  IERR
       REAL X, T, C, TAL, POP, WORK, QW, AVE
       LOGICAL ROBUST
       DIMENSION X(PLIM,N), T(P), CAT(P), C(PLIM,KMAX), 
     $  TAL(CATTOT,KMAX), Z(N), POP(KMAX), 
     $  IND(N), WORK(N), QW(N), AVE(N)
       
************************************************************************
*         Constants.  
************************************************************************
       INTEGER MAXIT, NP
       REAL VIS0, VIS9
       DOUBLE PRECISION PI, TOLPAR
       LOGICAL dump
       PARAMETER (
     $            dump = .true.,              ! turn on ticker tape
     $            MAXIT = 500,                 ! max iterations
     $            NP = 18,                     ! length of primes array
     $            PI = 3.141592653589793D0,    ! factor in Student pdf
     $            TOLPAR = 1D-4,               ! set convergence tolerance
     $            VIS0 = 0.001,                ! initial viscosity
     $            VIS9 = 0.500)                ! final viscosity
       
*               local variables
       INTEGER 
     $         G, H, i, J, J1, J2, K0, L,    ! counters
     $         HCUM,                         ! offset into TAL
     $         IALL,                         ! allocation error #
     $         KP,                           ! which prime constant
     $         LIMCAT, LIMCON,               ! array bounds
     $         PCAT, PCON,                   ! # categorical and continuous
     $         PRIME(NP),                    ! for shuffling
     $         RANK                          ! rank of correlation matrix
     
       REAL 
     $      ETA,                     ! numerical viscosity
     $      SEPS, SLIT,              ! machine constants
     $      Y,                       ! temp scalar
     $      VISPOW                   ! viscosity update exponent
     
       DOUBLE PRECISION
     $                  CTERM,       ! term in Student pdf
     $                  D,           ! temp scalar
     $                  DEPS, DLIT,  ! double precision machine constants
     $                  LAD,         ! log(abs(det(R)))
     $                  NU,          ! degrees of freedom of Student dist
     $                  NUPP,        ! nu+p
!     $                  ODK,         ! 1/k
     $                  ODN,         ! 1/n
     $                  QT,          ! sum of Q over L
     $                  TOL,         ! convergence tolerance
     $                  TOT          ! a total
     
*              allocatable arrays
       REAL, DIMENSION (:,:), ALLOCATABLE :: Q, QOLD, W
       DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE :: DQ, DS1, DS2, DT,
     $                                                  DEV, S
       DOUBLE PRECISION, DIMENSION (:,:), ALLOCATABLE :: R, U, V
       
*              external functions
       INTEGER HAT                 ! draw random #
       REAL RMACHINE
       DOUBLE PRECISION DMACHINE !!,
       LOGICAL DSAFDIV             ! safe to divide function
 
       SAVE PRIME
       DATA PRIME / 907, 911, 919, 929, 937, 941, 947, 953, 967, 
     $  971, 977, 983, 991, 997, 1009, 1013, 1019, 1021 /
       
************************************************************************
*                      begin       
************************************************************************
       IERR = 0 
       SEPS = RMACHINE(1)
       SLIT = RMACHINE(2)
       DEPS = DMACHINE(1)
       DLIT = DMACHINE(2)
       KP = HAT(NP)
       
       IF (K .LE. 1 .OR. K .GE. N) THEN   ! validate
         IERR = 3
         RETURN
       END IF
       
*          count categories
       PCON = 0
       PCAT = 0
       DO J = 1, P
         IF (CAT(J) .EQ. 0) THEN
           PCON = PCON + 1
         ELSE
           PCAT = PCAT + 1
         END IF
       END DO
       LIMCON = MAX(PCON, 1)
       LIMCAT = MAX(PCAT, 1)
       
*          request workspace
       ALLOCATE (DEV(LIMCON), DS1(LIMCON), DS2(LIMCON), DT(LIMCAT), 
     $  Q(N,KMAX), QOLD(N,KMAX), W(N,KMAX), DQ(KMAX), S(LIMCON), 
     $  R(LIMCON,LIMCON), U(LIMCON,LIMCON), V(LIMCON,LIMCON), STAT=IALL)
       IF (IALL .NE. 0) THEN
         PRINT *, 'tclust: trouble getting memory'
         IERR = 7
         RETURN
       END IF
       
       NU = DBLE(MAX(1, M - 1))               ! small degree of freedom
!!       NU = DBLE(MAX(1, N - KMAX - PCON))  ! large degree of freedom
       NUPP = NU + PCON
       ODN = 1. / DBLE(N)
       CTERM = ALGAMA(NUPP / 2.D0) - ALGAMA(NU / 2.D0) -
     $             0.5D0 * PCON * LOG(PI * NU)
!       TOL = TOLPAR * DBLE(N)
       TOL = 3 * N * KMAX * SEPS
       PRINT *, 'using tol of ', TOL
       ETA = VIS0
       VISPOW = (LOG(VIS9) / LOG(VIS0)) **(1. / FLOAT(MAXIT - 1))

*        set Q from initial partition
       DO i = 1, N
         L = Z(i)
         DO J = 1, K
           IF (J .EQ. L) THEN
             Q(i,J) = 1.
           ELSE
             Q(i,J) = 0.
           END IF
         END DO ! next j
       END DO ! next i
       
*        Uniform initial W
       DO i = 1, N
         DO L = 1, K
           W(i,L) = 1.
         END DO
       END DO
       
************************************************************************
*          main loop 
************************************************************************
       DO G = 1, MAXIT
       
*             quit?
         IF (G > ITER) GO TO 90       
         
*            populations
         DO L = 1, K
           CALL ADD (Q(1,L), N, POP(L))
         END DO
         
************************************************************************
*             delete empty clusters
************************************************************************
         K0 = K
         DO L = K0, 1, -1
           IF (POP(L) < SEPS) THEN
             IF (K .LE. KMIN) THEN
               DEALLOCATE (DEV, DS1, DS2, DT, Q, QOLD, W,
     $             DQ, S, R, U, V, STAT=IALL)
               IERR = 1
               RETURN
             END IF  ! too few?
             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
               Q(i,L) = Q(i,K)
               QOLD(i,L) = QOLD(i,K)
               W(i,L) = W(i,K)
             END DO
             POP(L) = POP(K)
             POP(K) = 0.
             K = K - 1
           END IF  ! empty?
         END DO  ! next L
                   
************************************************************************
*          let each center be a weighted average
************************************************************************
         DO L = 1, K
           DO H = 1, CATTOT            ! clear tally
             TAL(H,L) = 0.
           END DO
           HCUM = 0
           DO J = 1, P
             IF (CAT(J) .EQ. 0) THEN      ! continous and cyclic
               DO i = 1, N
                 WORK(i) = X(J,i)         ! data
                 QW(i) = Q(i,L) * W(i,L)  ! weight
               END DO
         
               IF (T(J) > SLIT) THEN        ! modulo arithmetic
                 CALL WCYMEN (WORK, QW, IND, N, T(J), AVE, i)
                 IF (i .EQ. 0) THEN
                   PRINT *, 'WCYMEN returns null answer!'
                   PRINT *, 'work: ', WORK
                   PRINT *, 'qw: ', QW
                   PRINT *, 'n: ', N
                   PRINT *, 't(j): ', T(J)
                   STOP
                 END IF
                 i = MOD(PRIME(KP) * L, i) + 1
                 C(J,L) = AVE(i) 
               ELSE IF (ROBUST) THEN ! scalar data
                 CALL WMED (WORK, QW, IND, N, C(J,L))
               ELSE                             ! weighted mean
                 CALL WMEAN (WORK, QW, N, C(J,L))
               END IF   ! cyclic ?
               
             ELSE
               C(J,L) = 0.
               DO i = 1, N
                 H = HCUM + NINT(X(J,i))
                 TAL(H,L) = TAL(H,L) + Q(i,L)
               END DO
               HCUM = HCUM + CAT(J)
             END IF ! continuous or categorical ?
             
           END DO     ! next j
         END DO    ! next L

*                  skip ahead if all categorical variables
         LAD = 0.D0
         IF (PCON < 1) GO TO 50

************************************************************************
*          deviations of continuous variables
************************************************************************
         DO J = 1, PCON
           DEV(J) = 0.D0
         END DO
         
         DO i = 1, N
           DO L = 1, K
             CALL DTPFC (X(1,i), P, T, CAT, CATTOT, LIMCON, LIMCAT, 
     $                    C(1,L), TAL(1,L), POP(L), SLIT, DS1, DT)
             DO J = 1, PCON
               DEV(J) = DEV(J) + Q(i,L) * DS1(J) **2
             END DO   ! next dimension
           END DO  ! next center
         END DO  ! next object
         
         DO J = 1, PCON
           DEV(J) = SQRT(DEV(J) * ODN)
         END DO

************************************************************************
*             compute the pooled correlation matrix
************************************************************************
         DO J2 = 1, PCON                    ! zero
           DO J1 = J2, PCON
             R(J1,J2) = 0.D0
           END DO 
         END DO
         TOT = 0.

         DO L = 1, K
           DO i = 1, N
           
!!!!!!!!!!!!!!!!
                IF (ISNAN(Q(i,L))) STOP 'NaN in Q at pooled correl!'
!!!!!!!!!!!!!!!!!!!!!!!!!!!                

             CALL DTPFC (X(1,i), P, T, CAT, CATTOT, LIMCON, LIMCAT, 
     $                    C(1,L), TAL(1,L), POP(L), SLIT, DS1, DT)
     
*                  standardize 
             DO J = 1, PCON
               IF (DSAFDIV(DS1(J), DEV(J))) DS1(J) = DS1(J) / DEV(J)
               IF (ISNAN(DS1(J))) STOP 'NaN in DS1!'
             END DO
             D = DPROD(Q(i,L), W(i,L)) 

*               accumulate             
             DO J2 = 1, PCON
               DO J1 = J2, PCON
                 R(J1,J2) = R(J1,J2) + D * (DS1(J1) * DS1(J2))                 
               END DO    ! next J1
             END DO    ! next J2
             TOT = TOT + D  ! Tyler & Vardi's formula
           END DO    ! next center
         END DO    ! next point         
         IF (DSAFDIV(1.D0, TOT)) THEN
           D = 1. / TOT
         ELSE
           D = ODN
         END IF

         DO J2 = 1, PCON
           DO J1 = J2, PCON
             R(J1,J2) = R(J1,J2) * D
             R(J2,J1) = R(J1,J2)         ! fill upper triangle
             IF (ISNAN(R(J1,J2))) STOP 'NaN in R!'
           END DO 
         END DO

*             pad the diagonal
        DO J = 1, PCON
          R(J,J) = R(J,J) + DEPS * PCON
        END DO

!!!!!!!!!!!!!!!              dump correl matrix         
!          PRINT *, 'using correlation matrix: '
  22      FORMAT (G14.4, $)
  23      FORMAT ()
!          DO J1 = 1, PCON  
!            DO J2 = 1, PCON
!              PRINT 22, R(J1,J2)
!            END DO
!            PRINT 23
!          END DO
!!!!!!!!!!!!!!!!!!!!!!!!          

*          pseudo-invert 
         CALL DRDI (R, PCON, RANK, LAD, U, S, V, DS1, DS2, IERR)
         IF (RANK < 1) THEN
           IERR = 4
           DEALLOCATE (DEV, DS1, DS2, DT, Q, QOLD, W,
     $             DQ, S, R, U, V, STAT=IALL)
           RETURN
         END IF

         IF (IERR .NE. 0) THEN
           IERR = 5
           DEALLOCATE (DEV, DS1, DS2, DT, Q, QOLD, W,
     $             DQ, S, R, U, V, STAT=IALL)
           RETURN
         END IF
         
  50     DO L = 1, K    ! save old memberships
           DO i = 1, N
             QOLD(i,L) = Q(i,L)
           END DO
         END DO

************************************************************************
*           membership Q each object each cluster
************************************************************************
         DO i = 1, N
           DO L = 1, K
             CALL DTPFC (X(1,i), P, T, CAT, CATTOT, LIMCON, LIMCAT,
     $                   C(1,L), TAL(1,L), POP(L), SLIT, DS1, DT)
              
*              standardize
             DO J = 1, PCON
               IF (DSAFDIV(DS1(J), DEV(J))) DS1(J) = DS1(J) / DEV(J)
             END DO

*               multivariate t log probability (Mahalanobis distance)
             DO J = 1, PCON ; DS2(J) = 0.D0 ; END DO
             DO J2 = 1, PCON
               DO J1 = 1, PCON
                 DS2(J1) = DS2(J1) + R(J1,J2) * DS1(J2)
               END DO
             END DO 
             D = 0.D0
             DO J = 1, PCON
               D = D + DS1(J) * DS2(J)
             END DO
             
             IF (D < 0.D0) THEN
               PRINT *, 'NEG mahalanobis dist! ', D
               D = 0.D0
!!!!!!!!!!!!!!!              dump correl matrix         
         PRINT *, 'spectrum:  ', S
          PRINT *, 'inverse correlation matrix: '
          DO J1 = 1, PCON  
           DO J2 = 1, PCON
              PRINT 22, R(J1,J2)
            END DO
            PRINT 23
          END DO
!!!!!!!!!!!!!!!!!!!!!!!!          
             END IF
             
             W(i,L) = SNGL(NUPP / (NU + D))      ! update weight 
             TOT = NUPP * LOG(1.D0 + D / NU)     ! Student log prob
             TOT = CTERM - 0.5D0 * (LAD + TOT)

*                   add in categorical log probabilities
             DO J = 1, PCAT
               DT(J) = LOG(DT(J))
             END DO
             CALL DADD (DT, PCAT, D)
             TOT = TOT + D
             DQ(L) = TOT                  ! argument to exponent
           END DO  ! next L
           
*             reduce magnitudes of log probabilities for numerical purposes
           D = DQ(1)
           DO L = 1, K
             D = MAX(D, DQ(L))
           END DO
           DO L = 1, K
             DQ(L) = DQ(L) - D
           END DO
             
*                  multivariate probability density
!!           ODK = 1.D0 / DBLE(K)
           DO L = 1, K
             DQ(L) = EXP(DQ(L)) * DBLE(POP(L)) * ODN  ! multiply by prior prob.
!!             DQ(L) = EXP(DQ(L)) * ODK                  ! uniform prior prob.
           END DO     ! next L
           CALL DADD (DQ, K, QT)
           
************************************************************************
*                 make a new cluster?
************************************************************************
           IF (QT .LE. 0.D0) THEN
             PRINT *, 'tclust: zero probability object: ', i
             IERR = 9
             RETURN           
           ELSE IF (QT < DLIT) THEN
             IF (K .GE. KMAX) THEN
               PRINT *, 'tclust:  too many clusters to insert'
               IERR = 8
               RETURN
             END IF
             DO L = 1, K                 ! erase memberships in other clusters
               Q(I,L) = 0.
             END DO
             K = K + 1
             DO H = 1, N
               IF (H .EQ. I) THEN
                 Q(H,K) = 1.
                 QOLD(H,K) = 1.
               ELSE
                 Q(H,K) = 0.
                 QOLD(H,K) = 0.
               END IF
             END DO
             DO H = 1, N
               W(H,K) = 1.
             END DO
             DO J = 1, P
               C(J,K) = X(J,I)
             END DO
             DO H = 1, CATTOT             ! clear tally
               TAL(H,K) = 0.
             END DO
             HCUM = 0
             DO J = 1, P                  ! count the lone object
               IF (CAT(J) > 0) THEN
                 H = HCUM + NINT(X(J,I))
                 TAL(H,K) = 1.
                 HCUM = HCUM + CAT(J)
               END IF
             END DO
             POP(K) = 1.
      
************************************************************************
*            update memberships for this object
************************************************************************
           ELSE IF (DSAFDIV(1.D0, QT)) THEN
             QT = 1.D0 / QT
             DO L = 1, K                  ! normalize probability
               Q(i,L) = SNGL(DQ(L) * QT)
             END DO
           ELSE
             PRINT *, 'tclust: cant divide by qt !'
             DO L = 1, K
               Q(i,L) = 1. / FLOAT(K)
             END DO
           END IF
         END DO         ! next i

************************************************************************
*               check for convergence
************************************************************************
         TOT = 0.D0
         DO L = 1, K
           DO i = 1, N
             TOT = TOT + (Q(i,L) - QOLD(i,L)) **2
           END DO
         END DO         
         IF (TOT < TOL) GO TO 90
         
*                 apply viscosity
         DO L = 1, K
           DO i = 1, N
             Q(i,L) = ETA * QOLD(i,L) + (1. - ETA) * Q(i,L)
           END DO
         END DO

         IF (dump) THEN
           PRINT *, 'tclust: iter ', G, ' clusters: ', K,  
     $      ' rank: ', RANK, ' at convergence test, tot: ', TOT
         END IF
     
         ETA = ETA **VISPOW
       END DO  ! next g

       IERR = 2      ! too many iterations
       DEALLOCATE (DEV, DS1, DS2, DT, Q, QOLD, W,
     $            DQ, S, R, U, V, STAT=IALL)
       RETURN
       
************************************************************************
*             set Z to maximum membership
************************************************************************
  90     DO i = 1, N
           Y = 0.
           DO L = 1, K
             IF (Q(i,L) > Y) THEN
               Y = Q(i,L)
               Z(i) = L
             END IF
           END DO
         END DO
                
*           modes
       DO L = 1, K
         HCUM = 0
         DO J = 1, P
           IF (CAT(J) > 0) THEN
             Y = 0.
             DO H = 1, CAT(J) 
               IF (TAL(HCUM+H,L) > Y) THEN
                 Y = TAL(HCUM+H,L)
                 C(J,L) = FLOAT(H)
               END IF
             END DO
             HCUM = HCUM + CAT(J)
           END IF
         END DO  ! next variable
       END DO  ! next cluster
         
!!!!!!!!!!!!!!!! experiment, show DEV
       IF (dump) THEN
         PRINT *, 'tclust: importance weights are: '
         DO J = 1, PCON
           PRINT 22, 1. / DEV(J)
         END DO 
         PRINT 23
       END IF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       
        
       DEALLOCATE (DEV, DS1, DS2, DT, Q, QOLD, W,
     $             DQ, S, R, U, V, STAT=IALL)
       IF (IALL .NE. 0) THEN
         PRINT *, 'tclust: trouble releasing memory!'
         IERR = 6
       END IF
       
       RETURN
      END ! of tclust
