C @(#)knnif.for	19.1 (ES0-DMD) 02/25/03 13:27:32
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.IDENTIFICATION
C
C  Program KNNIF
C
C.AUTHOR
C
C  F. Murtagh, ST-ECF, Garching.             Version 1.0  16 June 1986
C  F. Murtagh                                Version 2.0  Oct. 1988
C                                            (New std. ifs.)
C
C.PURPOSE
C
C  Pass parameters to, and execute, KNN (K-NNs Discriminant Analysis).
C
C.INPUT/OUTPUT
C
C  P1 to P4 contain parameters; these are: input table name, cardinality
C  of group 1 of the training set; test set table name; and number of
C  nearest neighbours.
C
C.ALGORITHM
C
C . uses Midas Table interface routines;
C . input table is assumed to contain entries in single precision;
C . no select or null values;
C . all columns of the input table are sent to the KNN routine;
C . description of the KNN  routine is in the corresponding program;
C . the output produced is to the screen only;
C . storage is limited only by the Midas Table system, and by the
C   overall system limitations;
C
C.MODIFICATIONS
C
C   NEW STANDARD INTERFACES, F. MURTAGH, AUG. 88
C   New Table File Interface M.Peron     SEP91
C
C-----------------------------------------------------------------------
        PROGRAM KNNIF
C
        CHARACTER*6     FORM
        CHARACTER*60     TRAIN, TEST
        INTEGER         NACTV,ISTAT,N1,N,M,K,TID1,NSORTC
        INTEGER         NAC,NAR1,IPTR1,TID2,I,N3,M3,IPTR2
        INTEGER         LEN,NSEL,IWPTR1,IWPTR2,NAR,NAR3,NTOT
        INTEGER         INULL,IADD0,IADD1,IADD2,IERR,N2
        INTEGER         MADRID, KUN, KNUL, DTYPE
        INCLUDE         'MID_INCLUDE:TABLES.INC/NOLIST'
        COMMON          /VMR/MADRID(1)
        INCLUDE         'MID_INCLUDE:TABLED.INC/NOLIST'
C
C
C ...   Assuming this is our first action ...
C
        CALL STSPRO('KNNIF')
        DTYPE = D_R4_FORMAT
C
C ...   Get input parameters, passed as keywords with the command.
C
        CALL STKRDC('P1',1,1,60,NACTV,TRAIN,KUN,KNUL,ISTAT)
        CALL STKRDI('INPUTI',1,1,NACTV,N1,KUN,KNUL,ISTAT)
        CALL STKRDC('P3',1,1,60,NACTV,TEST,KUN,KNUL,ISTAT)
        CALL STKRDI('INPUTI',2,1,NACTV,K,KUN,KNUL,ISTAT)
C
C ...   Read input tables.
C
        CALL TBTOPN(TRAIN,16,TID1,ISTAT)
        CALL TBIGET(TID1,M,N,NSORTC,NAC,NAR1,ISTAT)
        CALL TBCMAP(TID1,0,IPTR1,ISTAT)
C
        CALL TBTOPN(TEST,16,TID2,ISTAT)
        CALL TBIGET(TID2,M3,N3,NSORTC,NAC,NAR1,ISTAT)
        CALL TBCMAP(TID2,0,IPTR2,ISTAT)
C
C ...   Some error checking on input.
C
        IF (N.LT.1.OR.M.LT.1.OR.M3.LT.1.OR.N3.LT.1) THEN
           CALL STTPUT(' Nos. of rows/columns are less than 1.',ISTAT)
           CALL STTPUT(' What sort of a table is this ??',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
        IF (M3.NE.M) THEN
           CALL STTPUT
     X     (' Dimension of training and test sets different.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
        DO I = 1, M
           CALL TBFGET(TID1,I,FORM,LEN,DTYPE,ISTAT)
           IF (DTYPE.NE.D_R4_FORMAT) THEN
              CALL STTPUT(' Illegal format:',ISTAT)
              CALL STTPUT(' Only R*4 column type allowed.',ISTAT)
              CALL STTPUT(' Aborting.',ISTAT)
              GOTO 9000
           ENDIF
        ENDDO
        DO I = 1, M3
           CALL TBFGET(TID2,I,FORM,LEN,DTYPE,ISTAT)
           IF (DTYPE.NE.D_R4_FORMAT) THEN
              CALL STTPUT(' Illegal format:',ISTAT)
              CALL STTPUT(' Only R*4 column type allowed.',ISTAT)
              CALL STTPUT(' Aborting.',ISTAT)
              GOTO 9000
           ENDIF
        ENDDO
C
        CALL CHSEL(TID1,N,NSEL)
        IF (NSEL.NE.N) THEN
           CALL STTPUT(' Not all rows are SELECTed. ',ISTAT)
           CALL STTPUT(' In current implementation, MUST select all.',
     X                                                         ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
        CALL CHSEL(TID2,N3,NSEL)
        IF (NSEL.NE.N3) THEN
           CALL STTPUT(' Not all rows are SELECTed. ',ISTAT)
           CALL STTPUT(' In current implementation, MUST select all.',
     X                                                         ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
        CALL TDMGET(4*M*N,IWPTR1,ISTAT)
        CALL TBCMAP(TID1,1,IPTR1,ISTAT)
        CALL MAPSM(MADRID(IPTR1),MADRID(IWPTR1),
     1              NAR,N,M)
        CALL TDMGET(4*M3*N3,IWPTR2,ISTAT)
        CALL TBCMAP(TID2,1,IPTR2,ISTAT)
        CALL MAPSM(MADRID(IPTR2),MADRID(IWPTR2),
     1              NAR3,N3,M3)
        NTOT = N*M
        CALL CHNULL(MADRID(IWPTR1),NTOT,INULL)
        IF (INULL.NE.0) THEN
           CALL STTPUT
     X     (' Null entries in the table are not allowed.',ISTAT)
           CALL STTPUT
     X     (' Use SELECT, and then construct another table.',
     X                                                    ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
        NTOT = N3*M
        CALL CHNULL(MADRID(IWPTR2),NTOT,INULL)
        IF (INULL.NE.0) THEN
           CALL STTPUT
     X     (' Null entries in the table are not allowed.',ISTAT)
           CALL STTPUT
     X     (' Use SELECT, and then construct another table.',
     X                                                   ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
C ...   Allocate storage.
C
        N2 = N - N1
        CALL GETSTOR(K,IADD0)
        CALL GETSTOR(K,IADD1)
        CALL GETSTOR(K,IADD2)
C
C ...   We now have all the storage we need; APPL will look after
C       a few further odds and ends - e.g. setting a few hard-wired
C       options, rearranging output, etc.  Then it will call the
C       application program.
C
        IERR = 0
        CALL APPL(MADRID(IWPTR1),N1,N2,N,MADRID(IWPTR2),N3,M,K,
     X  MADRID(IADD0),MADRID(IADD1),MADRID(IADD2),IERR)
        IF (IERR.NE.0) THEN
           CALL STTPUT(' IERR not 0 on return from KNN.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
C ...   Finish up: write output file to disk, and release memory space.
C
        CALL RELSTOR(K,IADD0)
        CALL RELSTOR(K,IADD1)
        CALL RELSTOR(K,IADD2)
C
        CALL TBTCLO(TID1,ISTAT)
        CALL TBTCLO(TID2,ISTAT)
C
 9000   CONTINUE
        CALL STSEPI
C
        END
C----------------------------------------------------------------------
        SUBROUTINE APPL(TRAIN,N1,N2,N,TEST,N3,M,K,DK,
     X                          KLIST,KPOP,IERR)
C
        INTEGER N,M,K,N1,N2,N3,IERR
        REAL*4  TRAIN(N,M), TEST(N3,M), DK(K), KPOP(K), KLIST(K)
C
        CALL KNN(TRAIN,N1,N2,N,TEST,N3,M,K,KLIST,DK,KPOP,IERR)
        IF (IERR.NE.0) RETURN
C
        RETURN
        END
C-------------------------------------------------------------------------
        SUBROUTINE CHNULL(X,NT,NULL)
C ...   Check if null values are present.
        REAL*4  X(NT), TRNULL
        INTEGER TINULL,I,NULL,NT
        DOUBLE PRECISION TDNULL
        CALL TBMNUL(TINULL,TRNULL,TDNULL)
        NULL = 0
        DO I = 1, NT
           IF (X(I).EQ.TRNULL) NULL = NULL + 1
           IF (NULL.GT.0) RETURN
        ENDDO
        RETURN
        END
C-------------------------------------------------------------------------
        SUBROUTINE TBL_CHGP(X,NT,IL,IH)
C ...   Check if unacceptable group membership values are present.
        INTEGER NT,IL,IH,I
        REAL*4  X(NT)
        IL = 0
        IH = 0
        DO I = 1, NT
           IF (X(I).LT.0.0) IL = IL + 1
           IF (X(I).GT.2.0) IH = IH + 1
           IF (IL.GT.0.OR.IH.GT.0) RETURN
        ENDDO
        RETURN
        END
C-------------------------------------------------------------------------
        SUBROUTINE GETSTOR(NVALS,IPTR)
C ...   Allocate storage space for NVALS real values.
        INTEGER NVALS,IPTR,ISTAT
        CALL TDMGET(4*NVALS,IPTR,ISTAT)
        RETURN
        END
C-------------------------------------------------------------------------
        SUBROUTINE RELSTOR(NVALS,IPTR)
C ...   Release storage space (already allocated) for NVALS real values.
        INTEGER NVALS,IPTR,ISTAT
        CALL TDMFRE(4*NVALS,IPTR,ISTAT)
        RETURN
        END
