C @(#)ldaif.for	17.1.1.1 (ES0-DMD) 01/25/02 17:16:49
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 LDAIF
C
C.AUTHOR
C
C  F. Murtagh, ST-ECF, Garching.             Version 1.0  14 June 1986
C  F. Murtagh                                Version 2.0  Oct. 1988
C                                            (New std. ifs.)
C
C.PURPOSE
C
C  Pass parameters to, and execute, LDA (Fisher's Linear Discriminant
C  Analysis).
C
C.INPUT/OUTPUT
C
C  P1 and P2 contain parameters; these are: input table name, and output
C  table name.
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 LDA routine;
C . description of the LDA  routine is in the corresponding program;
C . the final column of the table input is taken as the definition of
C   the group memberships (values of 1, 2 and 0=unassigned are acceptable);
C . the output produced consists of a table with one column, giving the
C   projections on the axis connecting the two group centres of gravity;
C   value 0 divides the groups.
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. 1988
C  New Table File System    M. Peron    SEP91
C-----------------------------------------------------------------------
        PROGRAM LDAIF
C
        CHARACTER*6     FORM
        CHARACTER*60     NAMEIN, NAMEOUT
        CHARACTER*16    UNIT, LABEL
        INTEGER         NACTV,ISTAT,TID1,NCOLIN,NROW
        INTEGER         NSORTC,NAC,NAR,I,LEN,NSEL,NTOT
        INTEGER         IPTRIN,IWPTRIN,IPTRGP,INULL
        INTEGER         ILOW,IHIGH,IACTV,NCOLOUT
        INTEGER         TID,NNCOL,NNROW,NSC1,NAROUT
        INTEGER         IPTROUT,IWPTROUT,M
        INTEGER         IADD0,IADD1,IADD2,IADD3,IADD4,IADD5
        INTEGER         IADD6,IADD7,IERR
        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('LDAIF')
        DTYPE = D_R4_FORMAT
C
C ...   Get table name as a keyword passed in the command line.
C
        CALL STKRDC('P1',1,1,60,NACTV,NAMEIN,KUN,KNUL,ISTAT)
C
C ...   Read input table.
C
        CALL TBTOPN(NAMEIN,16,TID1,ISTAT)
        CALL TBIGET(TID1,NCOLIN,NROW,NSORTC,NAC,NAR,ISTAT)
C
C ...   Some error checking on input.
C
        IF (NROW.LT.1.OR.NCOLIN.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
        DO I = 1, NCOLIN
           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

        CALL CHSEL(TID1,NROW,NSEL)
        IF (NSEL.NE.NROW) 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
        NSEL = NROW
C
C ...   Pointer IPTRIN will map the beginning of the data in NAMEIN ...
        CALL TDMGET(4*NROW*NCOLIN,IWPTRIN,ISTAT)
        CALL TBCMAP(TID1,1,IPTRIN,ISTAT)
C ...   ... and pointer IPTRGP will map the last column (gp. memberships).
        CALL MAPSM(MADRID(IPTRIN),MADRID(IWPTRIN),
     1              NAR,NROW,NCOLIN-1)
        CALL TBCMAP(TID1,NCOLIN,IPTRGP,ISTAT)
        NTOT = NROW*(NCOLIN-1)
        CALL CHNULL(MADRID(IWPTRIN),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
        CALL TBL_CHGP(MADRID(IPTRGP),NROW,ILOW,IHIGH)
        IF (ILOW.NE.0.OR.IHIGH.NE.0) THEN
           CALL STTPUT(' Group memberships are unacceptable.',ISTAT)
           CALL STTPUT(' Check last column of the input table.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF

C
C
C ...   Get the output table name.
C
        CALL STKRDC('P2',1,1,60,IACTV,NAMEOUT,KUN,KNUL,ISTAT)
C
C ...   Now create the output (projections) table.
C
        NCOLOUT = 1
        CALL TBTINI(NAMEOUT,F_TRANS,17,NCOLOUT,NROW,TID,ISTAT)
        FORM = 'G14.6'
        UNIT = ' '
        LABEL = 'PROJNS'
        CALL TBCINI(TID,DTYPE,1,FORM,UNIT,LABEL,NSORTC,ISTAT)
C
        CALL TBIGET(TID,NNCOL,NNROW,NSC1,NAC,NAROUT,ISTAT)
        CALL TDMGET(4*NCOLOUT*NCOLIN,IWPTROUT,ISTAT)
        CALL TBCMAP(TID,1,IPTROUT,ISTAT)
        CALL MAPSM(MADRID(IPTROUT),MADRID(IWPTROUT),
     1              NAROUT,NROW,NCOLOUT)

C
C ...   Allocate storage.
C
        M = NCOLIN-1
        CALL GETSTOR(M,IADD0)
        CALL GETSTOR(2*M,IADD1)
        CALL GETSTOR(M*M,IADD2)
        CALL GETSTOR(M,IADD3)
        CALL GETSTOR(M,IADD4)
        CALL GETSTOR(M,IADD5)
        CALL GETSTOR(NROW*M,IADD6)
        CALL GETSTOR(NROW,IADD7)
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(NROW,M,MADRID(IWPTRIN),MADRID(IPTRGP),
     X  MADRID(IWPTROUT),MADRID(IADD0),MADRID(IADD1),MADRID(IADD2),
     X  MADRID(IADD3),MADRID(IADD4),
     X  MADRID(IADD5),MADRID(IADD6),MADRID(IADD7),IERR)
        CALL MAPBG(MADRID(IWPTROUT),MADRID(IPTROUT),
     1             NAROUT,NROW,NCOLOUT)
        IF (IERR.EQ.1) THEN
           CALL STTPUT(' IERR = 1 on return from LDA.',ISTAT)
           CALL STTPUT(' Group memberships are not 1 and 2.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
        IF (IERR.EQ.2) THEN
           CALL STTPUT(' IERR = 2 on return from LDA.',ISTAT)
           CALL STTPUT
     X(' Singular total varariance-covariance matrix.',ISTAT)
           CALL STTPUT
     X(' Is nrows > ncols? or is some column identically valued?',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
        IF (IERR.NE.0) THEN
           CALL STTPUT(' IERR not 0 on return from LDA.',ISTAT)
           CALL STTPUT(' Aborting.',ISTAT)
           GOTO 9000
        ENDIF
C
C ...   Finish up: write output file to disk, and release memory space.
C
        CALL RELSTOR(M,IADD0)
        CALL RELSTOR(2*M,IADD1)
        CALL RELSTOR(M*M,IADD2)
        CALL RELSTOR(M,IADD3)
        CALL RELSTOR(M,IADD4)
        CALL RELSTOR(M,IADD5)
        CALL RELSTOR(NROW*M,IADD6)
        CALL RELSTOR(M,IADD7)
C
        CALL TBTCLO(TID1,ISTAT)
C
        CALL TBIPUT(TID,NCOLOUT,NROW,ISTAT)
        CALL TBSINI(TID,ISTAT)
        CALL TBTCLO(TID,ISTAT)
C
 9000   CONTINUE
        CALL STSEPI
C
        END
C----------------------------------------------------------------------
        SUBROUTINE APPL(N,M,AIN,AGP,AOUT,MEAN,GPM,TOT,W1,W2,W3,
     X                  ACREA,GP,IERR)
C
C ...   This is now the application routine - note the setting of
C       parameters, which in a future epoch may be passed if desired
C       on the command line.
C
        INTEGER N,M,I,J,IERR,IPRINT
        REAL*4  AIN(N,M), AGP(N), AOUT(N), MEAN(M), GPM(2,M), TOT(M,M),
     X          W1(M), W2(M), W3(M), ACREA(N,M)
        INTEGER*4 GP(N)
C
        IPRINT = 3
C
        DO I = 1, N
           GP(I) = AGP(I)
           DO J = 1, M
              ACREA(I,J) = AIN(I,J)
           ENDDO
        ENDDO
C
        CALL LDA(N,M,ACREA,GP,IPRINT,MEAN,GPM,TOT,W1,W2,W3,IERR)
C
        IF (IERR.NE.0) RETURN
C
        DO I = 1, N
           AOUT(I) = ACREA(I,1)
        ENDDO
C
        RETURN
        END
C-------------------------------------------------------------------------
        SUBROUTINE CHNULL(X,NT,NULL)
C ...   Check if null values are present.
        REAL*4  X(NT), TRNULL
        INTEGER TINULL,NT,NULL,I
        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,I,IL,IH
        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