!C   ARGUMENTS    X      - INPUT MATRIX OF DIMENSION NV BY THE SUM OF 
!C                           ND(I) CONTAINING THE GROUP DATA.        
!C                           (SEE REMARK 1).                       
      REAL, ALLOCATABLE :: X(:,:)
      Integer :: IX, NG, NV
      Integer, ALLOCATABLE :: ND(:)

!C                IX     - INPUT ROW DIMENSION OF MATRIX X EXACTLY AS
!C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE
!C                           CALLING PROGRAM.                          
!C                NG     - INPUT NUMBER OF DISTINCT GROUPS.           
!C                NV     - INPUT NUMBER OF DATA VARIABLES.           
!C                ND     - INPUT VECTOR OF LENGTH NG.  ND(I) CONTAINS
!C                           THE NUMBER OF DATA VECTORS FROM THE I-TH
!C                           GROUP IN THE DATA MATRIX X.

      REAL, ALLOCATABLE :: P(:), XM(:,:)
      Integer :: IXM, IW, ID
      REAL, ALLOCATABLE :: W(:,:), D(:,:)
            
!C                P      - INPUT/OUTPUT VECTOR OF LENGTH NG. ON INPUT,      
!C                           P(I) CONTAINS THE BAYESIAN PRIOR PROBABILITY  
!C                           THAT A RANDOM OBSERVATION BELONGS TO THE     
!C                           I-TH GROUP. (SEE REMARK 2).                 
!C                         ON OUTPUT, P CONTAINS THE CONSTANTS OF THE   
!C                           NG LINEAR DISCRIMINANT FUNCTIONS.         
!C                XM     - OUTPUT MATRIX OF DIMENSION NV BY NG.       
!C                           THE I-TH COLUMN OF XM CONTAINS THE      
!C                           MEAN VECTOR FOR THE I-TH GROUP.        
!C                IXM    - INPUT ROW DIMENSION OF MATRIX XM EXACTLY AS         
!C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE      
!C                           CALLING PROGRAM.                                

!C                W      - OUTPUT MATRIX OF DIMENSION NV BY NG.  THE       
!C                           I-TH COLUMN OF W CONTAINS THE VECTOR FOR     
!C                           THE I-TH LINEAR DISCRIMINANT FUNCTION.      
!C                IW     - INPUT ROW DIMENSION OF MATRIX W EXACTLY AS   
!C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE    
!C                           CALLING PROGRAM.                              
!C                D      - OUTPUT PERFORMANCE MATRIX OF DIMENSION NG BY   
!C                           NG. (SEE REMARK 3).                         
!C                ID     - INPUT ROW DIMENSION OF MATRIX D EXACTLY AS   
!C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE 
!C                           CALLING PROGRAM.                           

!C   REMARKS  1.  ON INPUT, THE MATRIX X CONTAINS THE ND(1)+ND(2)+... 
!C                ...+ND(NG) SAMPLE (COLUMN) VECTORS. THE FIRST ND(1)
!C                COLUMNS OF MATRIX X SHOULD CONTAIN THE DATA FOR GROUP    
!C                ONE. THE NEXT ND(2) COLUMNS OF MATRIX X SHOULD CONTAIN  
!C                THE DATA FOR GROUP TWO, AND SO ON.                     

      REAL, ALLOCATABLE :: Y(:), Y1(:), Xtmp(:,:), S(:), Thres(:)   ! Working arrays
      Integer :: NThres
      Real :: Xtc

!      OPEN (1, FILE="data_test_disc.txt") !"Disc_Inp/JMA_VES_TC_All_1-5.txt")
      OPEN (1, FILE="Disc_Inp/JMA_VES_TC_All_1-5.txt")
      READ (1,*) NG, NV, N
      NThres=NG-1

      ALLOCATE (X(NV,N), ND(NG), P(NG), XM(NV,NG), W(NV,NG), D(NG,NG))
      ALLOCATE (Thres(NThres))
      ALLOCATE (Y(NV), Y1(N), Xtmp(N,NV), S(((NV+1)*NV)/2))

      Read (1,*) (Thres(i),i=1,NThres)

      print*,NG, NV, N, Thres(:)


      ND(:)=0
      Do i=1,N
         read (1,*) Xtc, (Xtmp(i,j),j=1,NV)
         Y1(i) = Group (Xtc)
         print "(20F7.3)",Y1(i),(Xtmp(i,j),j=1,NV)
         Do k=1,NG; if (Y1(i)==k) ND(k)=ND(k)+1 ; Enddo
      Enddo
      print*, (ND(k),k=1,NG)
      Close (1)
      ! Sort Y1:
      Do i=1, N-1
         Do j=i+1,N
            if (Y1(j)<Y1(i)) then
               Xtc=Y1(j); Y1(j)=Y1(i); Y1(i)=Xtc
               Y(:)=Xtmp(j,:); Xtmp(j,:)=Xtmp(i,:); Xtmp(i,:)=Y(:)
            endif
         Enddo
      Enddo

      ! Create Data array:
      Do iv=1, NV
         X(iv,:)=Xtmp(:,iv)
      Enddo

      ! PRIOR PROBABILITY:
      Do i=1,NG; P(i)=real(ND(i))/real(N); Enddo
      !Do i=1,NG; P(i)=1.0/Real(NG); Enddo

      ! Parameter:
      IX=NV; IXM=NV; IW=NV; ID=NG;  
      CALL ODNORM (X,IX,NG,NV,ND,P,XM,IXM,S,W,IW,D,ID,IER)

      ! Discriminal Equations:
      Do i=1, NG
         print*,"Group ",i
         print "(20F10.3)", P(i), (W(j,i),j=1,NV)
      Enddo

      ! Probability:
      Do i=1,NG
         print "(50F10.3)", (D(i,j),j=1,NG)
      Enddo

      ! Testing using Discr. functions:
      Do i=1,N
         Do k=1,NG
            Xtmp(i,k)=0
            Do j=1,NV
               Xtmp(i,k)=Xtmp(i,k)+W(j,k)*X(j,i) 
            Enddo
            Xtmp(i,k)=Xtmp(i,k)+P(k)
         Enddo
         Y(i)=GRoup_ID(Xtmp(i,1:NG),NG)
      Enddo
      
      ! Print test:
      i1=0
      Do k=1,NG
         Do j=1,ND(k)
            i1=i1+1
            print*,k, Int(Y(i1))
         Enddo
      Enddo

      CONTAINS
      ! -----------------------------------------------------------------
      Integer Function GRoup_ID (X,N)
         Integer :: N, i, j, ix
         Real :: X(N), Xmax
         Xmax=X(1); ix=1
         Do i=2,N
            if (X(i)>Xmax) ix=i
         Enddo
         GRoup_ID=ix
      End Function GRoup_ID
      ! -----------------------------------------------------------------
      Integer Function Group (x)
       ! Identify which group x is belong to
       Real :: x
       Integer :: k
       Do k=1,NG
          if (k==1) then
             if (x <= Thres(k)) Group=1
          elseif (k==NG) then
             if (x > Thres(k-1)) Group=NG
          else
             if (x > Thres(k-1).AND. x <= Thres(k)) Group=k
          endif
       Enddo
      End Function Group

      END  ! PROGRAM
     ! -----------------------------------------------------------------
     Include "FDA.for"