C @(#)knn.for	19.1 (ES0-DMD) 02/25/03 13:27:14
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 K-NN DISCRIMINANT ANALYSIS,                        C
C  with two groups defined in a training set, and with            C
C  assignment of members of a test set.                           C
C                                                                 C
C  Parameters:                                                    C
C                                                                 C
C  TRAIN(N,M)     training set, where first N1 rows relate to the C
C                 first group, and the next N2 rows to the second C
C                 group.  Must have N1 + N2 = N.                  C
C  TEST(N3,M)     test set;                                       C
C  K              number of nearest neighbours to consider;       C
C  KLIST(K), DK(K), KPOP(K)   are used for storing the K NNs,     C
C                 their distances to the object under consider-   C
C                 ation, and the group to which the NNs belong.   C
C                                                                 C
C  NOTE: this routine is written in VAX-11 Fortran and not in     C
C  Fortran-77.                                                    C
C  It also contains MIDAS STTPUT calls.                       C
C                                                                 C
C  F. Murtagh, ESA/ESO/STECF, Garching.  February 1986.           C
C                                                                 C
C-----------------------------------------------------------------C
        SUBROUTINE KNN(TRAIN,N1,N2,N,TEST,N3,M,K,KLIST,DK,KPOP,IERR)
        REAL TRAIN(N,M),TEST(N3,M),DK(K),D
        INTEGER KLIST(K),KPOP(K)
        INTEGER I,IX,N1,N2,N3,N,M,K,IERR,ISTAT
        INTEGER IND,I2,IK,IK2,NUM1,NUM2
        CHARACTER*80 LINE
C
        WRITE(LINE,20)
        CALL STTPUT(LINE,ISTAT)
   20   FORMAT(' Object  -->  group with probability')
        WRITE(LINE,30)
        CALL STTPUT(LINE,ISTAT)
   30   FORMAT(' ')
C
        DO I = 1, N3
        DO IX = 1, K
        KLIST(IX) = 0
        DK(IX) = 1.E+15
        KPOP(IX) = 0
        ENDDO
        IND = 0
                DO I2 = 1, N1
                CALL DISTK(I2,TRAIN,I,TEST,N,N3,M,D)
                IF (IND.GT.0) THEN
                   DO IK = 1, K
                        IF (D.LT.DK(IK)) THEN
                           IF (IK.LT.K) THEN
                              DO IK2 = K, IK+1, -1
                                 DK(IK2) = DK(IK2-1)
                                 KLIST(IK2) = KLIST(IK2-1)
                                 KPOP(IK2) = KPOP(IK2-1)
                              ENDDO
                           ENDIF
                           DK(IK) = D
                           KLIST(IK) = I2
                           KPOP(IK) = 1
                           GOTO 100
                        ENDIF
                    ENDDO
                ELSE
                    IND = IND + 1
                    DK(IND) = D
                    KLIST(IND) = I2
                    KPOP(IND) = 1
                ENDIF
  100           CONTINUE
                ENDDO
C
C
C
                DO I2 = N1+1, N
                CALL DISTK(I2,TRAIN,I,TEST,N,N3,M,D)
                IF (IND.GT.0) THEN
                   DO IK = 1, K
                        IF (D.LT.DK(IK)) THEN
                           IF (IK.LT.K) THEN
                              DO IK2 = K, IK+1, -1
                                 DK(IK2) = DK(IK2-1)
                                 KLIST(IK2) = KLIST(IK2-1)
                                 KPOP(IK2) = KPOP(IK2-1)
                              ENDDO
                           ENDIF
                           DK(IK) = D
                           KLIST(IK) = I2
                           KPOP(IK) = 2
                           GOTO 200
                        ENDIF
                    ENDDO
                ELSE
                    IND = IND + 1
                    DK(IND) = D
                    KLIST(IND) = I2
                    KPOP(IND) = 2
                ENDIF
  200           CONTINUE
                ENDDO
C
        NUM1 = 0
        NUM2 = 0
        DO IX = 1, K
           IF (KPOP(IX).EQ.1) NUM1 = NUM1 + 1
           IF (KPOP(IX).EQ.2) NUM2 = NUM2 + 1
        ENDDO
C       (Error check:)
        IF ((NUM1+NUM2).NE.K) THEN
           LINE(1:) = ' Total assignments to gp. 1 and to gp. 2 '
     X                //'dont equal K; is K too large ? '
           CALL STTPUT(LINE,ISTAT)
           IERR = 1
           RETURN
        ENDIF
        IF (NUM1.GT.NUM2) THEN
           WRITE (LINE,500) I,FLOAT(NUM1)*100./FLOAT(K)
           CALL STTPUT(LINE,ISTAT)
        ENDIF
        IF (NUM2.GT.NUM1) THEN
           WRITE (LINE,525) I,FLOAT(NUM2)*100./FLOAT(K)
           CALL STTPUT(LINE,ISTAT)
        ENDIF
        IF (NUM1.EQ.NUM2) THEN
           WRITE (LINE,550) I,FLOAT(NUM1)*100./FLOAT(K)
           CALL STTPUT(LINE,ISTAT)
        ENDIF
  500   FORMAT(I6,6X,'    1   ',F8.2,'%')
  525   FORMAT(I6,6X,'    2   ',F8.2,'%')
  550   FORMAT(I6,6X,' 1 or 2 ',F8.2,'%   (equiprobable).')
        ENDDO
C
        RETURN
        END
C
        SUBROUTINE DISTK(I,ARR1,J,ARR2,NR1,NR2,NC,D)
        INTEGER NC,NR1,NR2,K,J,I
        REAL ARR1(NR1,NC),ARR2(NR2,NC),D
C
        D = 0.0
        DO K = 1, NC
        D = D + (ARR1(I,K)-ARR2(J,K))**2
        ENDDO
C
        RETURN
        END
