C @(#)lda.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 Carry out a LINEAR DISCRIMINANT ANALYSIS, assigning ungrouped items C C to the closest group centre, using the Mahalanobis distance. C C C C C C To call: CALL LDA(N,M,DATA,GP,IPRINT,MEAN,MGP,TOTAL,DIFFVEC, C C W1,W2,IERR) where C C C C C C N, M : integer dimensions of ... C C DATA : input data (real). C C On output the first column of DATA contains projections of C C all N items on the line connecting the two group means. C C Zero is the boundary point. C C GP : Integer vector of length N giving group assignments. An C C unassigned item has group 0. Otherwise groups 1 and 2 will C C be looked for, and other values here are not acceptable. C C IPRINT: integer; print options (= 3: full; otherwise none). C C MEAN : real vector of length M (number of attributes or variables). C C MGP : real array of dimensions 2 by M. C C TOTAL : real array of dimensions M by M; on output contains inverse C C of total variance/covariance matrix. C C DIFFVEC: real vector of length M. C C W1, W2: real vectors of length M. C C C C C C Inputs here are N, M, DATA, GP, IPRINT (and IERR). C C The principle output information is contained in DATA. C C IERR = 1 means that more than two groups have been specified; IERR C C = 2 means that the total variance-covariance matrx was singular. C C C C Note: we require N > M > 2, to prevent the seeking of the inverse C C of a singular matrix. C C C C C C F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986. C C C C-----------------------------------------------------------------------C SUBROUTINE LDA(N,M,DATA,GP,IPRINT,MEAN,MGP,TOTAL,DIFFVEC,W1,W2, X IERR) REAL DATA(N,M), TOTAL(M,M), MEAN(M), MGP(2,M), DIFFVEC(M) REAL W1(M), W2(M),D INTEGER GP(N), NOG(2), G,IMAT INTEGER N,M,IPRINT,IERR,NEFF,I,J INTEGER K,J1,J2 C C Form global mean. C IERR = 0 NEFF = 0 DO 50 I = 1, N IF (GP(I).NE.0) NEFF = NEFF + 1 IF (GP(I).LE.2) GOTO 40 IERR = 1 GOTO 9000 40 CONTINUE 50 CONTINUE C DO 200 J = 1, M MEAN(J) = 0.0 DO 100 I = 1, N IF (GP(I).NE.0) MEAN(J) = MEAN(J) + DATA(I,J) 100 CONTINUE MEAN(J) = MEAN(J)/FLOAT(NEFF) 200 CONTINUE C C Form (total) variance-covariance matrix. C DO 500 J1 = 1, M DO 400 J2 = 1, M TOTAL(J1,J2) = 0.0 DO 300 I = 1, N IF (GP(I).NE.0) TOTAL(J1,J2) = TOTAL(J1,J2) + X (DATA(I,J1)-MEAN(J1))*(DATA(I,J2)-MEAN(J2)) 300 CONTINUE TOTAL(J1,J2) = TOTAL(J1,J2)/FLOAT(NEFF) 400 CONTINUE 500 CONTINUE C IMAT = 1 C IF (IPRINT.EQ.3) CALL OUTMTL(IMAT,M,TOTAL) C C Form group means. C DO 700 J = 1, M DO 600 K = 1, 2 MGP(K,J) = 0.0 600 CONTINUE 700 CONTINUE C DO 900 I = 1, N G = GP(I) IF (G.EQ.0) GOTO 900 NOG(G) = NOG(G) + 1 DO 800 J = 1, M MGP(G,J) = MGP(G,J) + DATA(I,J) 800 CONTINUE 900 CONTINUE C DO 1100 K = 1, 2 DO 1000 J = 1, M MGP(K,J) = MGP(K,J)/NOG(K) 1000 CONTINUE 1100 CONTINUE C C Invert variance-covariance matrix. C CALL MINVL(M,TOTAL,D,W1,W2) IF (ABS(D).GT.0.00001) GOTO 1150 IERR = 2 GOTO 9000 1150 CONTINUE C C Form difference vector of group mean vectors. C DO 1200 J = 1, M DIFFVEC(J) = MGP(1,J) - MGP(2,J) MEAN(J) = (MGP(1,J) + MGP(2,J))/2. 1200 CONTINUE C C Determine projections and output them. C CALL PROJXL(N,M,DATA,MEAN,W1,TOTAL,DIFFVEC) IF (IPRINT.EQ.3) CALL OUTPXL(N,M,DATA) C C C 9000 CONTINUE RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Output a matrix. C C C C-----------------------------------------------------------------------C SUBROUTINE OUTMTL(IMAT,M,ARRAY) C C Output array. C INTEGER IMAT,M,K1,K2 REAL ARRAY(M,M) C C IF (IMAT.EQ.1)WRITE (6,900) DO 100 K1 = 1, M WRITE (6,1000) (ARRAY(K1,K2),K2=1,M) 100 CONTINUE C C 900 FORMAT(' VARIANCE/COVARIANCE MATRIX FOLLOWS.',/) 1000 FORMAT(5(2X,F8.4)) RETURN END C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Invert a symmetric matrix and calculate its determinant. C C C C C C To call: CALL MATINV(M,ARRAY,DET,W1,W2) where C C C C C C M : dimension of ... C C ARRAY : input matrix which is replaced by its inverse. C C NORDER : degree of matrix (order of determinant) C C DET : determinant of input matrix. C C W1, W2 : work vectors of dimension M. C C C C C C Reference: Philip B Bevington, "Data Reduction and Error Analysis C C for the Physical Sciences", McGraw-Hill, New York, 1969, C C pp. 300-303. C C C C C C F. Murtagh, ST-ECF, Garching-bei-Muenchen, January 1986. C C C C----------------------------------------------------------------------C SUBROUTINE MINVL(M,ARRAY,DET,IK,JK) INTEGER M,K,I,J,L REAL ARRAY(M,M), IK(M), JK(M) REAL DET,SAVE,AMAX C DET = 1.0 DO 100 K = 1, M C Find largest element ARRAY(I,J) in rest of matrix. AMAX = 0.0 21 DO 30 I = K, M DO 30 J = K, M IF (ABS(AMAX)-ABS(ARRAY(I,J))) 24,24,30 24 AMAX = ARRAY(I,J) IK(K) = I JK(K) = J 30 CONTINUE C Interchange rows and columns to put AMAX in ARRAY(K,K). IF (AMAX) 41,32,41 32 DET = 0.0 GOTO 140 41 I = IK(K) IF (I-K) 21,51,43 43 DO 50 J = 1, M SAVE = ARRAY(K,J) ARRAY(K,J) = ARRAY(I,J) 50 ARRAY(I,J) = -SAVE 51 J = JK(K) IF (J-K) 21,61,53 53 DO 60 I = 1, M SAVE = ARRAY(I,K) ARRAY(I,K) = ARRAY(I,J) 60 ARRAY(I,J) = -SAVE C Accumulate elements of inverse matrix. 61 DO 70 I = 1, M IF (I-K) 63,70,63 63 ARRAY(I,K) = -ARRAY(I,K)/AMAX 70 CONTINUE DO 80 I = 1, M DO 80 J = 1, M IF (I-K) 74,80,74 74 IF (J-K) 75,80,75 75 ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K)*ARRAY(K,J) 80 CONTINUE DO 90 J = 1, M IF (J-K) 83,90,83 83 ARRAY(K,J) = ARRAY(K,J)/AMAX 90 CONTINUE ARRAY(K,K) = 1.0/AMAX 100 DET = DET * AMAX C Restore ordering of matrix. DO 130 L = 1, M K = M - L + 1 J = IK(K) IF (J-K) 111,111,105 105 DO 110 I = 1, M SAVE = ARRAY(I,K) ARRAY(I,K) = -ARRAY(I,J) 110 ARRAY(I,J) = SAVE 111 I = JK(K) IF (I-K) 130,130,113 113 DO 120 J = 1, M SAVE = ARRAY(K,J) ARRAY(K,J) = -ARRAY(I,J) 120 ARRAY(I,J) = SAVE 130 CONTINUE 140 RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Output projections of row points. C C C C-----------------------------------------------------------------------C SUBROUTINE OUTPXL(N,M,PRJN) INTEGER N,M,NUM,K,J REAL PRJN(N,M) C C NUM = 1 WRITE (6,1000) WRITE (6,1010) WRITE (6,1020) DO 100 K = 1, N WRITE (6,1030) K,(PRJN(K,J),J=1,NUM) 100 CONTINUE C 1000 FORMAT(1X,'PROJECTIONS OF ROW-POINTS FOLLOW.',/) 1010 FORMAT(' OBJECT PROJN') 1020 FORMAT(' ------ ------') 1030 FORMAT(I5,2X,F8.4) RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Form projections of row-points on factors. C C C C-----------------------------------------------------------------------C SUBROUTINE PROJXL(N,M,DATA,MEAN,VEC,TOTINV,DIFF) INTEGER N,M,K,L,I,J1,NUM,J2 REAL DATA(N,M),MEAN(M),VEC(M),TOTINV(M,M),DIFF(M) C NUM = 1 DO 300 K = 1, N DO 50 L = 1, M VEC(L) = DATA(K,L) 50 CONTINUE DO 200 I = 1, NUM DATA(K,I) = 0.0 DO 100 J1 = 1, M DO 75 J2 = 1, M DATA(K,I) = DATA(K,I) + (VEC(J1)-MEAN(J1))* X TOTINV(J1,J2)*DIFF(J2) 75 CONTINUE 100 CONTINUE 200 CONTINUE 300 CONTINUE C RETURN END