! Partitioning ! DIMENSION A(18,16), MEMGP(18), NUMGP(18), COMP(2) ! DIMENSION GPCEN(2,16) Integer, Parameter :: N=36, M=7, NG=5 DIMENSION A(N,M), MEMGP(N), NUMGP(N), COMP(NG) DIMENSION GPCEN(NG,M) OPEN(UNIT=21,STATUS='OLD',FILE='spectr.dat') !N = 18 !M = 16 !NG = 2 DO I = 1, N READ(21,1000) (A(I,J),J=1,M) 1000 FORMAT(8F7.1) ENDDO IERR = 0 NGP0 = 1 NTRIES = 15 ISEED = 37519 DO I = 1, NTRIES CALL RANDP(N,NG,MEMGP,ISEED) ! In the following, we optimise either by assigning (using ! minimum distances to group centres) or by exchanging. CALL MINDST(A,N,M,MEMGP,NGP0,NUMGP,GPCEN,NG,COMP,CTOT,ITER,IERR) !CALL EXCH(A,N,M,MEMGP,NGP0,NUMGP,GPCEN,NG,COMP,CTOT,ITER,IERR) IF (IERR.EQ.2) THEN WRITE(6,*) ' A GROUP HAS LESS THAN',NGP0,' MEMBERS.', & ' REDUCE THE NUMBER OF GROUPS AND TRY AGAIN.' STOP ENDIF IF (IERR.EQ.1) THEN WRITE(6,*) ' AN INVALID GROUP ASSIGNMENT HAS BEEN DETECTED', & ' (LESS THAN 1 OR GREATER THAN THE NO. OF GPS.).', & ' IS THE NUMBER OF GROUPS CORRECTLY SPECIFIED?' STOP ENDIF IF (IERR.NE.0) THEN WRITE(6,*) ' IERR =',IERR STOP ENDIF WRITE(6,*) ' SUM OF VARIANCES:',CTOT WRITE(6,*) ' NO. OF ITERATIONS:',ITER WRITE(6,*) ' ASSIGNMENTS:', (MEMGP(K),K=1,N) ENDDO END !--------------------------------------------------------------------- ! C ! Generate a random partition. C ! C ! Parameters: C ! C ! MEMGP(N) Group memberships, C ! NG number of groups, C ! ISEED seed for random number generator. C ! C ! Note: random number generator is machine-dependent. C ! C ! F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C ! C !-----------------------------------------------------------------C SUBROUTINE RANDP(N,NG,MEMGP,ISEED) DIMENSION MEMGP(N) DO 100 I = 1, N MEMGP(I) = 1 100 CONTINUE IF (NG.LE.1.OR.N.LE.1) GOTO 500 X = 1.0/FLOAT(NG) DO 400 I = 1, N VAL = RAN(ISEED) BNDRY = X ICL = 1 200 IF (ICL.EQ.NG) GOTO 300 IF (VAL.LT.BNDRY) GOTO 300 BNDRY = BNDRY + X ICL = ICL + 1 GOTO 200 300 MEMGP(I) = ICL 400 CONTINUE 500 CONTINUE RETURN END !--------------------------------------------------------------------- ! C ! Optimise the variances of a set of groups, by assigning C ! the objects in groups such that they are minimally distant C ! from group centres. C ! C ! Parameters: C ! C ! N, M, NG Numbers of rows, columns, groups, C ! A(N,M) initial data, C ! MEMGP(N) group memberships, C ! NGP0 minimum acceptable group cardinality, C ! NUMGP(NG) cardinalities of groups, C ! GPCEN(NG,M) group centres, C ! COMP(NG) compactness values for the groups, C ! CTOT sum of these compactnesses, C ! IERR error indicator (should be zero). C ! C ! IERR = 1: invalid group number (<1 or >NG), - is number of C ! groups correctly specified? IERR = 2: a group has < minimum C ! allowed number of members, - reduce the number of groups and C ! try again. C ! C ! F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C ! C !-----------------------------------------------------------------C SUBROUTINE MINDST(A,N,M,MEMGP,NGP0,NUMGP,GPCEN,NG, COMP,CTOT,ITER,IERR) DIMENSION A(N,M), MEMGP(N), NUMGP(NG), GPCEN(NG,M), COMP(NG) BIG = 1.0E+30 ONE = 0.999 CMAX = BIG ITER = 0 100 ITER = ITER + 1 IF (ITER.GT.15) GOTO 500 CALL GMEANS(A,N,M,MEMGP,NGP0,NUMGP,GPCEN,NG,IERR) CALL COMPCT(A,N,M,NG,MEMGP,GPCEN,COMP,CTOT) IF (IERR.NE.0) GOTO 500 IF (NG.LE.1) GOTO 500 IF (CTOT.GE.CMAX) GOTO 500 CMAX = CTOT*ONE DO 400 I = 1, N X = BIG DO 300 K = 1, NG Y = 0.0 DO 200 J = 1, M DIFF = GPCEN(K,J) - A(I,J) Y = Y + DIFF*DIFF 200 CONTINUE IF (Y.GE.X) GOTO 300 X = Y ICL = K 300 CONTINUE MEMGP(I) = ICL 400 CONTINUE GOTO 100 500 RETURN END !--------------------------------------------------------------------- ! C ! Optimise the variances of a set of groups, by exchanging C ! the objects between groups such that they are minimally C ! distant from group centres. C ! C ! Parameters: C ! C ! N, M, NG Numbers of rows, columns, groups, C ! A(N,M) initial data, C ! MEMGP(N) group memberships, C ! NGP0 minimum acceptable group cardinality, C ! NUMGP(NG) cardinalities of groups, C ! GPCEN(NG,M) group centres, C ! COMP(NG) compactness values for the groups, C ! CTOT sum of these compactnesses, C ! IERR error indicator (should be zero). C ! C ! IERR = 1: invalid group number (<1 or >NG), - is number of C ! groups correctly specified? IERR = 2: a group has < minimum C ! allowed number of members, - reduce the number of groups and C ! try again. C ! C ! F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C ! C !-----------------------------------------------------------------C !--------------------------------------------------------------------- SUBROUTINE EXCH(A,N,M,MEMGP,NGP0,NUMGP,GPCEN,NG, COMP,CTOT,ITER,IERR) DIMENSION A(N,M), MEMGP(N), NUMGP(NG), GPCEN(NG,M), COMP(NG) BIG = 1.0E+3 ONE = 0.999 CALL GMEANS(A,N,M,MEMGP,NGP0,NUMGP,GPCEN,NG,IERR) CALL COMPCT(A,N,M,NG,MEMGP,GPCEN,COMP,CTOT) IF (IERR.NE.0) GOTO 800 IF (NG.LE.1) GOTO 800 ITER = 0 I = 0 IS = 0 100 IS = IS + 1 IF (IS.GT.N) GOTO 800 200 I = I + 1 IF (I.LE.N) GOTO 300 ITER = ITER + 1 IF (ITER.GT.15) GOTO 800 I = 1 300 ICL = MEMGP(I) NUM = NUMGP(ICL) IF (NUM.LE.NGP0) GOTO 100 V = NUM EQ = BIG DO 600 K = 1, NG X = 0.0 DO 400 J = 1, M DIFF = GPCEN(K,J) - A(I,J) X = X + DIFF*DIFF 400 CONTINUE IF (K.NE.ICL) GOTO 500 FRAC1 = V/(V-1.0) EP = X*FRAC1 GOTO 600 500 FRAC2 = NUMGP(K) FRAC = FRAC2/(FRAC2+1.0) EK = FRAC*X IF (EK.GE.EQ) GOTO 600 EQ = EK IQ = K W = FRAC2 600 CONTINUE IF (EQ.GE.EP*ONE) GOTO 100 IS = 0 COMP(ICL) = COMP(ICL) - EP COMP(IQ) = COMP(IQ) + EQ CTOT = CTOT - EP + EQ FRAC1 = 1.0/(V-1.0) FRAC2 = 1.0/(W+1.0) DO 700 J = 1, M VAL = A(I,J) GPCEN(ICL,J) = (V*GPCEN(ICL,J)-VAL)*FRAC1 GPCEN(IQ,J) = (W*GPCEN(IQ,J)+VAL)*FRAC2 700 CONTINUE MEMGP(I) = IQ NUMGP(ICL) = NUM - 1 NUMGP(IQ) = NUMGP(IQ) + 1 GOTO 200 800 CONTINUE RETURN END !--------------------------------------------------------------------- ! C ! Standardize to zero mean and unit standard deviation. C ! C ! Parameters: C ! C ! N, M, NG Numbers of rows, columns, groups, C ! A(N,M) initial data, replaced by standardized values.C ! C ! F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C ! C !-----------------------------------------------------------------C SUBROUTINE STND(A,N,M) DIMENSION A(N,M) DO 500 J = 1, M X = 0.0 DO 100 I = 1, N X = X + A(I,J) 100 CONTINUE XBAR = X/FLOAT(N) X = 0.0 DO 200 I = 1, N DIFF = A(I,J) - XBAR X = X + DIFF*DIFF 200 CONTINUE IF (X.LE.0.0) X = 1.0 X = 1.0/SQRT(X) DO 300 I = 1, N A(I,J) = X*(A(I,J)-XBAR) 300 CONTINUE 500 CONTINUE RETURN END !--------------------------------------------------------------------- ! C ! Determine means of the groups. C ! C ! Parameters: C ! C ! N, M, NG Numbers of rows, columns, groups, C ! A(N,M) initial data, C ! MEMGP(N) group memberships, C ! NGP0 minimum acceptable group cardinality, C ! NUMGP(NG) cardinalities of groups, C ! GPCEN(NG,M) group centres, C ! IERR error indicator (should be zero). C ! C ! F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C ! C !-----------------------------------------------------------------C SUBROUTINE GMEANS(A,N,M,MEMGP,NGP0,NUMGP,GPCEN,NG,IERR) DIMENSION A(N,M), MEMGP(N), NUMGP(NG), GPCEN(NG,M) DO 200 K = 1, NG NUMGP(K) = 0 DO 100 J = 1, M GPCEN(K,J) = 0.0 100 CONTINUE 200 CONTINUE DO 500 I = 1, N ICL = MEMGP(I) IF (ICL.GE.1.AND.ICL.LE.NG) GOTO 300 IERR = 1 RETURN 300 CONTINUE NUMGP(ICL) = NUMGP(ICL) + 1 DO 400 J = 1, M GPCEN(ICL,J) = GPCEN(ICL,J) + A(I,J) 400 CONTINUE 500 CONTINUE DO 800 K = 1, NG NUM = NUMGP(K) IF (NUM.GE.NGP0) GOTO 600 IERR = 2 RETURN 600 CONTINUE X = 1.0/FLOAT(NUM) DO 700 J = 1, M GPCEN(K,J) = GPCEN(K,J)*X 700 CONTINUE 800 CONTINUE RETURN END !--------------------------------------------------------------------- ! C ! Determine compactness of the groups (i.e. their variances). C ! C ! Parameters: C ! C ! N, M, NG Numbers of rows, columns, groups, C ! A(N,M) initial data, C ! MEMGP(N) group memberships, C ! GPCEN(NG,M) group centres, C ! COMP(NG) variances of groups (output), C ! CTOT sum of these variances. C ! C ! F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C ! C !-----------------------------------------------------------------C SUBROUTINE COMPCT(A,N,M,NG,MEMGP,GPCEN,COMP,CTOT) DIMENSION A(N,M), MEMGP(N), GPCEN(NG,M), COMP(NG) CTOT = 0.0 DO 100 K = 1, NG COMP(K) = 0.0 100 CONTINUE DO 300 I = 1, N ICL = MEMGP(I) X = 0.0 DO 200 J = 1, M DIFF = GPCEN(ICL,J) - A(I,J) X = X + DIFF*DIFF 200 CONTINUE COMP(ICL) = COMP(ICL) + X CTOT = CTOT + X 300 CONTINUE RETURN END !---------------------------------------------------------------------