C @(#)part.for 17.1.1.1 (ES0-DMD) 01/25/02 17:16:38 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Generate a random partition. C C C C Parameters: C C C C MEMGP(N) Group memberships, C C NG number of groups, C C ISEED seed for random number generator. C C C C Note: random number generator is machine-dependent. C C C C F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C C C C-----------------------------------------------------------------C SUBROUTINE RANDP(N,NG,MEMGP,ISEED) INTEGER N,NG,ISEED,I,ICL INTEGER MEMGP(N) REAL X,VAL,BNDRY,SRAN C DO 100 I = 1, N MEMGP(I) = 1 100 CONTINUE C IF (NG.LE.1.OR.N.LE.1) GOTO 500 X = 1.0/FLOAT(NG) DO 400 I = 1, N VAL = SRAN(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 C 500 CONTINUE RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Uniform random number generator, Spaeth p. 143 C C Returns real in interval [0,1) C C C C-----------------------------------------------------------------C REAL FUNCTION SRAN(IY) INTEGER IY,IA,IC,M2,M,MIC,ITWD DOUBLE PRECISION HALFM, DATAN, DSQRT REAL S DATA M2/0/,ITWD/2/ IF (M2.NE.0) GOTO 20 M=1 10 M2=M M=ITWD*M2 IF (M.GT.M2) GOTO 10 HALFM=M2 IA=8*IDINT(HALFM*DATAN(1.D0)/8.D0) + 5 IC=2*IDINT(HALFM*(0.5D0-DSQRT(3.D0)/6.D0)) + 1 MIC=(M2-IC)+M2 S=0.5/HALFM 20 IY=IY*IA IF(IY.GT.MIC) IY=(IY-M2)-M2 IY=IY+IC IF (IY/2.GT.M2) IY=(IY-M2)-M2 IF(IY.LT.0) IY=(IY+M2)+M2 SRAN=FLOAT(IY)*S RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Optimise the variances of a set of groups, by assigning C C the objects in groups such that they are minimally distant C C from group centres. C C C C Parameters: C C C C N, M, NG Numbers of rows, columns, groups, C C A(N,M) initial data, C C MEMGP(N) group memberships, C C NGP0 minimum acceptable group cardinality, C C NUMGP(NG) cardinalities of groups, C C GPCEN(NG,M) group centres, C C COMP(NG) compactness values for the groups, C C CTOT sum of these compactnesses, C C IERR error indicator (should be zero). C C C C IERR = 1: invalid group number (<1 or >NG), - is number of C C groups correctly specified? IERR = 2: a group has < minimum C C allowed number of members, - reduce the number of groups and C C try again. C C C C F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C C C C-----------------------------------------------------------------C SUBROUTINE MINDST(A,N,M,MEMGP,NGP0,NUMGP,GPCEN,NG, X COMP,CTOT,ITER,IERR) REAL A(N,M), GPCEN(NG,M),COMP(NG),BIG,ONE REAL CMAX,X,Y,DIFF,CTOT INTEGER MEMGP(N),NUMGP(N),ITER,I,J,K,ICL INTEGER N,M,IERR,NG,NGP0 C 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 C 500 RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Optimise the variances of a set of groups, by exchanging C C the objects between groups such that they are minimally C C distant from group centres. C C C C Parameters: C C C C N, M, NG Numbers of rows, columns, groups, C C A(N,M) initial data, C C MEMGP(N) group memberships, C C NGP0 minimum acceptable group cardinality, C C NUMGP(NG) cardinalities of groups, C C GPCEN(NG,M) group centres, C C COMP(NG) compactness values for the groups, C C CTOT sum of these compactnesses, C C IERR error indicator (should be zero). C C C C IERR = 1: invalid group number (<1 or >NG), - is number of C C groups correctly specified? IERR = 2: a group has < minimum C C allowed number of members, - reduce the number of groups and C C try again. C C C C F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C C C C-----------------------------------------------------------------C SUBROUTINE EXCH(A,N,M,MEMGP,NGP0,NUMGP,GPCEN,NG, X COMP,CTOT,ITER,IERR) REAL A(N,M),GPCEN(NG,M), COMP(NG) REAL BIG,ONE,V,EQ,X,FRAC1,FRAC2,EP REAL EK,DIFF,FRAC,VAL,CTOT,W INTEGER MEMGP(N),NUMGP(NG),N,M,IERR,NG INTEGER ITER,I,ICL,NUM,IS,J,K,NGP0,IQ C BIG = 1.0E+30 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 C 800 CONTINUE RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Standardize to zero mean and unit standard deviation. C C C C Parameters: C C C C N, M, NG Numbers of rows, columns, groups, C C A(N,M) initial data, replaced by standardized values.C C C C F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C C C C-----------------------------------------------------------------C SUBROUTINE STND(A,N,M) INTEGER I,J,M,N REAL A(N,M),DIFF,XBAR,X C 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 C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Determine means of the groups. C C C C Parameters: C C C C N, M, NG Numbers of rows, columns, groups, C C A(N,M) initial data, C C MEMGP(N) group memberships, C C NGP0 minimum acceptable group cardinality, C C NUMGP(NG) cardinalities of groups, C C GPCEN(NG,M) group centres, C C IERR error indicator (should be zero). C C C C F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C C C C-----------------------------------------------------------------C SUBROUTINE GMEANS(A,N,M,MEMGP,NGP0,NUMGP,GPCEN,NG,IERR) REAL A(N,M), GPCEN(NG,M),X INTEGER MEMGP(N),NUMGP(NG),K,NG,I,M,N INTEGER ICL,J,NUM,IERR,NGP0 C DO 200 K = 1, NG NUMGP(K) = 0 DO 100 J = 1, M GPCEN(K,J) = 0.0 100 CONTINUE 200 CONTINUE C 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 C 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 C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Determine compactness of the groups (i.e. their variances). C C C C Parameters: C C C C N, M, NG Numbers of rows, columns, groups, C C A(N,M) initial data, C C MEMGP(N) group memberships, C C GPCEN(NG,M) group centres, C C COMP(NG) variances of groups (output), C C CTOT sum of these variances. C C C C F. Murtagh, ESA/ESO/STECF, Garching-bei-Muenchen, Feb. 1986. C C C C-----------------------------------------------------------------C SUBROUTINE COMPCT(A,N,M,NG,MEMGP,GPCEN,COMP,CTOT) REAL A(N,M), GPCEN(NG,M), COMP(NG),X REAL DIFF,CTOT INTEGER MEMGP(N) INTEGER N,M,NG,I,J,K,ICL C CTOT = 0.0 DO 100 K = 1, NG COMP(K) = 0.0 100 CONTINUE C 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 C RETURN END