C++*********************************************************************
C
C SUBKMEANS.F
C
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2005  Health Research Inc.                      *
C=*                                                                    *
C=* HEALTH RESEARCH INCORPORATED (HRI),                                *   
C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455.                   *
C=*                                                                    *
C=* Email:  spider@wadsworth.org                                       *
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 the *
C=* 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 GNU  *
C=* General Public License for more details.                           *
C=*                                                                    *
C=* You should have received a copy of the GNU General Public License  *
C=* along with this program; if not, write to the                      *
C=* Free Software Foundation, Inc.,                                    *
C=* 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.      *
C=*                                                                    *
C **********************************************************************
C
C SUBKMEANS
C
C	W():	WEIGHTS TO BE USED PER OBJECT.
C	INUM():	FACTORS TO BE USED TO BE USED.
C	NFAC:	NUMBER OF FACTORS TO BE USED
C	NFTOT:	TOTAL NUMBER OF FACTORS USED FOR MULTIVARIATE STATISTICAL
C	 	ANALYSIS (MSA) I.E; CORAN OR PCA.
C	NOBJ	NUMBER OF OBJECTS.
C       KCLASS   : NUMBER OF CLUSTERS
C       X(M,L)   : COORDINATES OF EACH OBJECT
C       IC(M)    : CLUSTER TO WHICH OBJECT M BELONGS.
C       CM(L,N)  : COORDINATES OF THE CENTROIDS OF EACH CLUSTER.
C       WSS(N)   : SUM OF THE SQUARED DISTANCES OF ALL OBJECTS
C                  (IN A CLUSTER) TO THEIR CENTROID.
C       NC(N)    : NUMBER OF OBJECTS IN EACH CLUSTER.
C       D        : SUM OF E(N)
C                  (THIS HAS TO BE REDUCED AFTER EACH NEW CLUSTERING)
C
C	CALLS:
C	  NEWKMEANS, FILSEQP, FILGET, SAVDN1, SAVD, SAVDC
C
C	NOTE: FILE ON UNIT=LUNF IS ALREADY OPEN.
C             CLOSE THE FILE ON RETURN FROM THIS ROUTINE.
C
C--*********************************************************************

        SUBROUTINE SUBKMEANS(X, CM, NIMI,IC, NC, WSS, AT, NOBJ, NFAC,
     &                       KCLASS, W, NFTOT, COO, INUM, BB, WW, DB,
     &                       IX_DIM,USE_DISK,ISEED,LUNF,LUND,ITYPE)

        INCLUDE 'CMBLOCK.INC' 
        INCLUDE 'CMLIMIT.INC' 

        DIMENSION    X(NFAC,IX_DIM), CM(NFAC,KCLASS)
        DIMENSION    IC(NOBJ), NIMI(NOBJ), NC(KCLASS), WSS(KCLASS)
        DIMENSION    AT(NFAC),W(NFAC), COO(NFTOT), INUM(NFAC)

        PARAMETER    (NLIST = 3)
        DIMENSION    DLIST(NLIST)

        LOGICAL               :: USE_DISK
        CHARACTER(LEN=MAXNAM) :: FINPAT,FINP
        INTEGER               :: ILIST(1)

C	SQUARE ROOTS OF THE WEIGHTS.
        DO J = 1, NFAC
           W(J) = SQRT(W(J))
        ENDDO

	IF (ISEED .EQ. 0)  THEN
C	   IC(NOBJ): ASSIGN A CLUSTER TO EACH OBJECT. 
           I = 1
           DO WHILE (I .LE. NOBJ)
              K = 1
              DO WHILE (K .LE. KCLASS)
                 IC(I) = K
                 K     = K + 1
                 I     = I + 1
              ENDDO
           ENDDO
	ELSE
C          RANDOM ASSIGNMENT
	   T = 1.0/FLOAT(KCLASS)
	   I = 1
	   DO WHILE (I .LE. NOBJ)
	      Q  = RAND_P(ISEED)
	      QT = T
	      K  = 1
	      DO WHILE (QT .LE. Q)
	         QT = QT +T 
	         K  = K + 1
	      ENDDO
	      IC(I) = K
	      I     = I + 1
	   ENDDO
	ENDIF

        IDR = 0
        CALL NEWKMEANS(NOBJ,NFAC,X,IC,KCLASS,CM,WSS,NC,D,IDR,
     &               INUM,NIMI,W,COO,NFTOT,IX_DIM,USE_DISK,LUNF,ITYPE)

C       CALCULATE CITERIA ...

c$omp   parallel do if(nfac.ge.64), private(i)
	DO I = 1,NFAC
           AT(I) = 0.0
        ENDDO

        DO L = 1,NOBJ
c$omp     parallel do if(nfac.ge.64), private(i)
          DO I = 1,NFAC
            AT(I) = AT(I) + X(I,L)
          ENDDO
        ENDDO

c$omp   parallel do if(nfac.ge.64), private(i)
        DO I = 1,NFAC
          AT(I) = AT(I)/NOBJ
        ENDDO

        BB = 0.0
c$omp   parallel do if(nfac.ge.64), private(i,j),reduction(+:bb)
        DO I = 1,NFAC
           DO J = 1,KCLASS
              BB = BB + (CM(I,J) - AT(I))**2
           ENDDO
        ENDDO

C	REWIND THE UNIT=LUNF

        IF (USE_DISK) THEN
           IF (ITYPE .EQ. 1) THEN
              CALL REW(LUNF,1)
           ELSE
              CALL REWF(LUNF,1)
           ENDIF
        ENDIF

	DO I=1,KCLASS
 	   WSS(I) = 0.0
	ENDDO

        DO L = 1, NOBJ
          IF (USE_DISK) THEN
             IF (ITYPE .EQ. 1) THEN
                READ(LUNF)   (COO(J), J = 1, NFTOT)
             ELSE
                READ(LUNF,*) (COO(J), J = 1, NFTOT)
             ENDIF

             DO J = 1, NFAC
                X(J,IX_DIM) = W(J) * COO(INUM(J))
             ENDDO
             K = IX_DIM
          ELSE
             K = L
          ENDIF

          WW = 0.0
c$omp     parallel do if(nfac.ge.64), private(i),reduction(+:ww)
          DO I = 1, NFAC
              WW = WW + (X(I,K) - CM(I,IC(L)))**2
          ENDDO
	  WSS(IC(L)) = WSS(IC(L))+WW
        ENDDO
        WW = 0.0
	DO I=1,KCLASS
 	   WW = WW + WSS(I)
	ENDDO

C       Davis-Bouldin criterium
	DB = 0.0
	DO J=1,KCLASS
	   RK = -1.0
	   DO I=1,KCLASS
	      IF (I.NE.J)  THEN
	         SMK = 0.0
c$omp            parallel do if(nfac.ge.64), private(l),reduction(+:smk)
	         DO L=1,NFAC
 	            SMK = SMK + (CM(L,J) - CM(L,I))**2
	         ENDDO
	         RJK = (WSS(J) + WSS(I) ) / SMK
	         RK  = AMAX1(RK,RJK)
	      ENDIF
	   ENDDO
 	   DB = DB + RK
	ENDDO

	DB = DB / K

C	WRITE THE OBJECTS IN EACH CLASS TO A FILE.
        NMAX = 0
        CALL FILSEQP(FINPAT,NLET,ILIST,NMAX,NIMA,
     &                'TEMPLATE FILENAME (e.g.: SEL***)',IRTFLG)

C	LOOKING THROUGH IC(), PUT ALL OBJECTS BELONGING TO THE SAME
C	CLASS IN A FILE. 

        DO ICLASS = 1, KCLASS

C	  CREATE A FILENAME FOR CLASS ICLASS
          CALL FILGET(FINPAT,FINP,NLET,ICLASS,IRTFLG)

C	  NUMBER OF OBJECT IN THIS GROUP
          IGRP = 0
          NLS  = 2
	  DO I = 1, NOBJ
            IF (ICLASS .EQ. IC(I)) THEN

C             THE OBJECT BELONG TO CLASS ICLASS, STORE IT.
              IGRP     = IGRP + 1
              DLIST(1) = IGRP
              DLIST(2) = NIMI(I)
              IAP      = 0
              CALL SAVDN1(LUND,FINP,DLIST,NLS,IGRP-1,IAP)
            ENDIF
	  ENDDO

C         CLOSE EACH CLASS FILE OPENED BY SAVDN1 ROUTINE.
          CLOSE(LUND)
	ENDDO

        DO I = 1, NOBJ
           DLIST(1) = I
           DLIST(2) = NIMI(I)
           DLIST(3) = IC(I)
           CALL SAVD(LUND, DLIST, NLIST, IRTFLG)
        ENDDO

        CALL  SAVDC
        CLOSE(LUND)
        CLOSE(LUNF)

        END