! Linear Discriminal Analysis

    PROGRAM LDA
      Implicit None
      INTEGER :: NGROUP, NVAR, NMAX, N, NThres
          ! So nhom, So bien, Dung luong MAX moi nhom, Cac nguong de chia nhom
      INTEGER, ALLOCATABLE :: NG(:)    ! Dung luong mau thuc te

      REAL, ALLOCATABLE :: X(:,:,:)     ! So lieu (OBS, VAR, GROUP)
      REAL, ALLOCATABLE :: Y(:), Y1(:), Xtmp(:,:)   ! Working arrays
      REAL, ALLOCATABLE :: MEAN(:,:), CMEAN(:), Thres(:) 
          ! Trung binh tung nhom va Trung binh chung
      REAL, ALLOCATABLE :: CG(:,:,:), C(:,:), C1(:,:)
          ! Covariance Matrixes
      REAL, ALLOCATABLE :: P(:), COEF(:,:), Tmp(:)
      REAL, ALLOCATABLE :: F(:,:), X1(:,:)
      INTEGER, ALLOCATABLE :: LOBS(:), LFSC(:)
      Character (Len=30) :: Fmt
      Real :: Xtc, Temp
      Integer :: i,j,k,k1,i1, L, j1, j2, N1
      ! Doc SL
      OPEN (1, FILE="data_test_disc.txt") !"Disc_Inp/JMA_VES_TC_All_1-5.txt")
      READ (1,*) NGROUP, NVAR, N              ! Num of groups, Num of variables and Total samples
      NThres=NGROUP-1

      ALLOCATE( NG(NGROUP), MEAN(NGROUP, NVAR), CMEAN(NVAR) ) 
      ALLOCATE( CG(NVAR, NVAR, NGROUP), C(NVAR, NVAR),C1(NVAR,NVAR))
      ALLOCATE( P(NGROUP), COEF(NGROUP,0:NVAR), Tmp(NVAR) )
      Allocate (Y(N), Y1(N), Xtmp(N,Nvar), Thres(NThres))

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

      ! Reading origin data:
      NG(:)=0
      Do i=1,N
         read (1,*) Xtc, (Xtmp(i,j),j=1,Nvar)
         Y1(i) = Group (Xtc)
         Do k=1,NGROUP; if (Y1(i)==k) NG(k)=NG(k)+1 ; Enddo
      Enddo
      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
      ! Max size of groups' sample size:
      NMAX  = NG(1)
      Do K=2,NGROUP    ! Dung luong mau CUC DAI cua cac nhom
         IF (NG(K) > NMAX) NMAX = NG(K)
      Enddo
      ! 3D-Data Array: 
      ALLOCATE( X(NMAX, NVAR, NGROUP) )

      !Do K=1,NGROUP
      !   k1=0
      !   Do i=1,N
      !      if (Y1(i)==K) then
      !         k1=k1+1; X(k1,:,K) = Xtmp(i,:)
      !         print*,K,K1,Y1(i),NG(k) 
      !      endif
      !   Enddo
      !Enddo

      ! Re-store 2D-data array into 3D-array:
      i1=0
      Do K=1,NGROUP
         Do i=1,NG(K)
            i1=i1+1; X(i,:,k)=Xtmp(i1,:)
            print*, i1,i
         Enddo
      Enddo

      Open (3, file="Disc_Inp/Test.txt",status="Unknown")

      write (3,"(50I5)") NGROUP, NVAR, (NG(k),k=1,NGROUP)
      Do K=1,NGROUP
         i=0
         Do j=1,NG(k)
            i=i+1; write (3,"(50F8.3)") (X(i,k1,K),k1=1,NVAR)
         Enddo
      Enddo
      Close (3)

      !N = 0   ! Tong dung luong mau
      !Do i=1,Ngroup
      !   N = N + NG(i)
      !Enddo

      ALLOCATE( F(N, NGROUP), X1(N, NVAR), LOBS(N), LFSC(N) )
      ! Xac dinh cac nhom ban dau
      L = 0
      Do K=1,NGROUP
         Do I=1,NG(K)
            L = L + 1
            LOBS(L) = K
         Enddo
      Enddo
      Do J = 1, NVAR
         CMEAN(J) = 0.  ! Trung binh chung
         Do k = 1, Ngroup
            MEAN(K,J) = 0.  ! Trung binh tung nhom
            Do i=1, NG(K)
               MEAN(K,J) = MEAN(K,J) + X(I,J,K)
            Enddo
            CMEAN(J) = CMEAN(J) + MEAN(K,J)
            MEAN(K,J) = MEAN(K,J)/REAL(NG(K))
         Enddo
         CMEAN(J) = CMEAN(J) /REAL(N)
      Enddo
      ! Tinh do lech
      Do K=1,Ngroup
         Do J=1, NVAR
            Do I=1,NG(K)
               X(I,J,K) = X(I,J,K) - CMEAN(J)
            Enddo
         Enddo
      Enddo
      Do J1=1,NVAR
         Do J2=1,NVAR
            C(J1,J2) = 0.
            Do K=1,Ngroup
               CG(J1,J2,K) = 0.
               Do I=1,NG(K)
                  CG(J1,J2,K) = CG(J1,J2,K) + & 
                        X(I,J1,K)*X(I,J2,K)
               Enddo
               CG(J1,J2,K) = CG(J1,J2,K)/REAL(NG(K))
               C(J1,J2) = C(J1,J2) + CG(J1,J2,K)*NG(K)
            Enddo
            C(J1,J2) = C(J1,J2)/REAL(N)
         Enddo
      Enddo

      ! Tinh xac suat tien nghiem
      Do K=1, Ngroup
         P(K) = REAL(NG(K))/REAL(N)
      Enddo

      ! Tinh ma tran nghich dao cua C
      C1(:,:) = C(:,:)
      CALL MINVERT(C1,NVAR)

      ! Tinh cac he so cua ham phan biet
      Do K=1, NGroup
         Do J1=1,NVAR
            COEF(K,J1) = 0.
            Do J2=1,NVAR
               COEF(K,J1) = COEF(K,J1) + MEAN(K,J2)*C1(J2,J1)
            Enddo
         Enddo
         COEF(K,0) = 0.
         DO J1=1,NVAR
            Tmp(J1) = 0.
            Do J2=1,NVAR
               Tmp(J1) = Tmp(J1) + MEAN(K,J2)*C1(J2,J1)
            Enddo
         Enddo
         Do J2=1,NVAR
            COEF(K,0) = COEF(K,0) + Tmp(J2)*MEAN(K,J2)
         Enddo

         COEF(K,0) = -0.5*COEF(K,0) + LOG(P(K))
         print*,'Cac he so ham phan biet ', K
         print '(8F9.4)',COEF(K,0:)
      Enddo

      ! Tinh ham phan biet cho tap mau phu thuoc
      Do J=1,NVAR
         N1 = 0
         Do K=1,Ngroup
            Do I=1,NG(K)
               N1 = N1+1
               X1(N1,J) = X(I,J,K) + CMEAN(J)
            Enddo
         Enddo
      Enddo
      Do I=1,N1
         Do K=1,NGROUP
            F(I,K) = 0.
            Do J=1,NVAR
               F(I,K) = F(I,K) + COEF(K,J)*X1(I,J)
            Enddo
            F(I,K) = F(I,K) + COEF(K,0)
         Enddo
         LFSC(I) = 1
         Temp = F(I,1)
         Do K=2,NGROUP
            IF (F(I,K)>Temp) THEN
               Temp = F(I,K)
               LFSC(I) = K
            EndIf
         Enddo
      Enddo
      write (Fmt,"(A,I2,A)") "(I5,",NGROUP,"F10.4,2I4)"
      Do I=1,N1
         print Fmt, I,(F(I,K),K=1,NGROUP),LOBS(I), LFSC(I)
      Enddo
    CONTAINS
    ! -----------------------------------------------------------------
    Integer Function Group (x)
      Real :: x
      Integer :: k
      Do k=1,NGROUP
         if (k==1) then
            if (x <= Thres(k)) Group=1
         elseif (k==NGROUP) then
            if (x > Thres(k-1)) Group=NGROUP
         else
            if (x > Thres(k-1).AND. x <= Thres(k)) Group=k
         endif
      Enddo
    End Function Group

    END  ! PROGRAM
    ! -----------------------------------------------------------------
    SUBROUTINE MINVERT(A,N)
      ! C/TRINH TINH MA TRAN NGHICH DAO CUA MA TRAN A
      ! INPUT: + A(N,N) MANG So lieu 
      !        + N SO HANG/COT CUA MA TRAN A
      ! OUTPUT: A MA TRAN NGHICH DAO CUA A
      !
      Integer :: N, M, i, j, k
      REAL :: A(N,N), R
      
      M=N-1
      DO I=1,N
         R=A(I,1)
         DO K=1,M
            A(I,K)=A(I,K+1)/R
         ENDDO
         A(I,N)=1.0/R
         DO J=1,N
            IF (I/=J) THEN
               R=A(J,1)
               DO K=1,M
                  A(J,K)=A(J,K+1)-R*A(I,K)
               ENDDO
               A(J,N)=-R*A(I,N)
            ENDIF
         ENDDO
      ENDDO 
      RETURN
   END SUBROUTINE 
   ! -----------------------------------------------------------------