PROGRAM t_logistic
      Use NetCDF
      USE Logistic_Regression
      IMPLICIT NONE

      Real, parameter    :: Miss=-99.0
      Integer, Parameter :: Nmon=12, Nbin=10
      Real, parameter    :: Bin(Nbin)=(/0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95/)
      Character (Len=255) :: Fname, Inp_Dir, Mask_File, Out_Dir
      Integer :: Ntime,M, i,j,k, ilead, ndf, ier, IO
      Real    :: aa,bb,cc,dd, scales
      !Real    :: PC, POD, TS, H1
      Real,Allocatable :: XX(:,:)
      Real (dp), Allocatable :: X(:,:), beta(:), se_beta(:)
      Integer, Allocatable   :: LG(:),N(:),NG(:),s1(:)
      Real (dp) :: tmp, chisq, devnce, tail_prob
      Logical :: Status
      !------------------------------------------------------------------------------
      Character (Len=10) :: Dr_type
      Real, Allocatable :: Regs(:,:), Lon(:), Lat(:)
      Real, Allocatable, Dimension (:,:,:,:) :: Dr_obs, Dr_fcs   ! (Nlon,Nlat,Nlead,Ntime) 
      Real, Allocatable, Dimension (:,:,:)   :: PC, POD, TS, HSS
      Real, Allocatable, Dimension (:,:,:,:) :: betas
      Real, Allocatable, Dimension (:,:,:,:) :: P0, P1           ! (Nlon,Nlat,Nlead,Nbin)
      Integer :: Nlon,Nlat,Nlead, ireg,ilon,ilat, Nt,it, icnt, icount(6), ic1(6), ii
      Integer :: Yr1,Yr2,Nmods,Ncase,Ngroup,ilen,idr
      Real    :: Dr_Thres
      Integer, Parameter :: NDr_lev=4, Nlen=3
      Character (Len=3)  :: Dr_lev(NDr_lev)=(/"mil","mod","sev","ext"/)
      !------------------------------------------------------------------------------
      Real    :: Px0(Nbin), Px1(Nbin), xtmp
      !------------------------------------------------------------------------------
      Namelist /File_Dir_Names/ Inp_Dir, Mask_File, Out_Dir !, TVar, RVar
      Namelist /Data_Info/      Yr1, Yr2, Dr_Thres, Dr_type, Nmods, Ncase
      !------------------------------------------------------------------------------
      Open (1, file="Namelist_Dr_Chars_Discr.txt", Status="Old")     ! Tinh voi APH 1961-2015
      Read (1, File_Dir_Names)
      Read (1, Data_Info)
      Close(1)
      Ntime=(Yr2-Yr1+1)*Nmon   ! Maximum timeseries
      Nlead=5                  ! All leadtimes
      Call rd_ncdf_2d (Mask_File,Lon,Lat,Nlon,Nlat,Regs,"reg") 

      Allocate (Dr_obs(Nlon,Nlat,Nlead,Ntime),Dr_fcs(Nlon,Nlat,Nlead,Ntime))
      Allocate (PC(Nlon,Nlat,Nlead),POD(Nlon,Nlat,Nlead),TS(Nlon,Nlat,Nlead),HSS(Nlon,Nlat,Nlead))
      Allocate (P0(Nlon,Nlat,Nlead,Nbin),P1(Nlon,Nlat,Nlead,Nbin))
      Dr_obs=Miss; Dr_fcs=Miss;
      PC=Miss; POD=Miss; TS=Miss; HSS=Miss;
      P0=Miss; P1=Miss

      ilen=1            ! Drought timescale: 1: 1-month, 3: 3-month, 6: 6-month
      idr =1            ! Drought level:     1: Moderate, 2: Severe, 3:Extreme
      !write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Dr_4_Logistic_Moderate_Len_",ilen,".dat"
      write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Dr_evt_and_model_out_4_discr_logis_Dr_"//trim(Dr_lev(idr))//"_len_",ilen,".dat"

      Open (1,file=Fname,status="Old",form='unformatted')

      IO=0; icnt=0; icount(:)=0; ic1(:)=0
      Do while (io==0)
         !if (IO/=0) exit
         Read (1,IOSTAT=io) ilead, ireg ! ilon, ilat
         if (IO/=0) exit
         Read (1) Nt, M, Ngroup
         icnt=icnt+1
         ii=mod(ilead,Nlead); if (ii==0) ii=Nlead
         icount(ii)=icount(ii)+1
print*,ilead,ireg,Nt, M, Ngroup; pause
         if (allocated (LG)) Deallocate (LG); if (allocated (X))  Deallocate (X)
         if (allocated (N))  Deallocate (N);  if (allocated (s1)) Deallocate (s1)
         if (allocated (beta)) Deallocate (beta)
         if (allocated (se_beta)) Deallocate (se_beta)
         if (allocated (XX)) Deallocate (XX)
         if (allocated (NG)) Deallocate (NG)
         Allocate (LG(Nt), X(Nt,M), N(Nt), s1(Nt))
         Allocate (beta(0:M),se_beta(0:M))
         Allocate (XX(Nt,M))
         Allocate (NG(Ngroup))
         if (icnt==1) Allocate (betas(Nlon,Nlat,Nlead,M+1))

         N(:) = 1
         Read (1) NG
         Read (1) LG
         Read (1) XX
         if (IO/=0) exit
!      Enddo
!print*,icount; stop
!print*,LG(1:Nt)
print*,XX(1:50,1)
         X(:,:)=XX(:,:); Deallocate (XX)
         CALL logistic(Nt,X,M,LG,N, chisq, devnce, ndf, beta, se_beta, ier)

         IF (ier /= 0) THEN
            ic1(ii)=ic1(ii)+1
            WRITE(*, *) 'Error number', ier
         ELSE
     !       WRITE(*, '(a, f9.3, a, f9.3, a, i4, a)')  &
     !         ' Deviance = ', devnce, '   Chi-squared = ', chisq,   &
     !         ' with ', ndf, ' deg. of freedom'
     !       tail_prob = 1.0_dp - chi_squared(ndf, chisq)
     !       IF (tail_prob < 0.01_dp) THEN
     !          WRITE(*, *) '*** Significantly bad fit at the 1% level ***'
     !       END IF
     !       WRITE(*, '(a, g12.4/)') ' Chi-squared tail prob. = ', tail_prob
     !       WRITE(*, *) '        Coefficient   Std.error'
     !       DO i = 0, M
     !          WRITE(*, '(" ", I3, "  ", 2g13.5)') i, beta(i), se_beta(i)
     !       END DO

            ! Rescale the beta: From t_lgstc3.f90 
            !---------------------------------------------------------------
            scales=1.0
            Do j=0,M
               scales=scales*ABS(beta(j))
            Enddo
            scales=SQRT(scales); beta = beta/scales
            !---------------------------------------------------------------
            aa=0; bb=0;cc=0;dd=0
            Px0(1:Nbin)=0; Px1(1:Nbin)=0
            Do i=1,Nt
               tmp=beta(0)
               Do j=1,M
                  tmp=tmp+beta(j)*x(i,j)         ! F = exp(b0 + b1.X1 + b2.X2 + ... + bk.Xk)
               Enddo
               tmp=EXP(tmp); tmp=tmp/(1+tmp);    ! p = F/(1 + F)
               if (tmp>0.030) then; s1(i)=1; else; s1(i)=0; endif
               Dr_obs(ilon,ilat,ilead,i)=LG(i)
               Dr_fcs(ilon,ilat,ilead,i)=tmp

               Do j=1,Nbin
                  !if (LG(i)==1.AND.tmp>=Bin(j)-0.05.AND.tmp<Bin(j)+0.05) Px1(j)=Px1(j)+1
                  !if (LG(i)==0.AND.tmp>=Bin(j)-0.05.AND.tmp<Bin(j)+0.05) Px0(j)=Px0(j)+1
                  if (tmp>=Bin(j)-0.05.AND.tmp<Bin(j)+0.05) then
                     if (LG(i)==1) Px1(j)=Px1(j)+1; if (LG(i)==0) Px0(j)=Px0(j)+1
                  endif
               Enddo

!               print"(2I4,F6.3)",LG(i),s1(i),tmp
               !if (LG(i)==s1(i)) print*,LG(i),s1(i),tmp
               if (s1(i)==1.AND.LG(i)==1) aa=aa+1
               if (s1(i)==1.AND.LG(i)==0) bb=bb+1
               if (s1(i)==0.AND.LG(i)==1) cc=cc+1
               if (s1(i)==0.AND.LG(i)==0) dd=dd+1
            Enddo
            Px1(:)=Px1(:)/Real(Nt); Px0(:)=Px0(:)/Real(Nt)
            P0(ilon,ilat,ilead,1:Nbin)=Px0(1:Nbin); P1(ilon,ilat,ilead,1:Nbin)=Px1(1:Nbin)
            print"(A3,10F8.3)","Bin",Bin(1:Nbin)
            print"(A3,10F8.3)","Px0",Px0(1:Nbin)
            print"(A3,10F8.3)","Px1",Px1(1:Nbin)
            print*,"   PC          POD "
            print"(4F8.3)",(aa+dd)/(aa+bb+cc+dd),aa/(aa+cc); 
            !pc=aa+bb+cc+dd; pc=(aa+dd)/pc
            PC(ilon,ilat,ilead) =(aa+dd)/(aa+bb+cc+dd)
            POD(ilon,ilat,ilead)=aa/(aa+cc); 
            TS(ilon,ilat,ilead) =aa/(aa+bb+cc); 
            HSS(ilon,ilat,ilead)=2.0*(aa*dd-bb*cc)/((aa+cc)*(cc+dd)+(aa+bb)*(bb+dd))
            !print*,"PC, POD, TS, HSS=",APC, POD, TS, AHSS
            !print*,POD, TS
         END IF
if (mod(icnt,30)==0) then; print*,"Lead=",ilead; pause; endif
      Enddo
print*,icount(:)
print*,ic1(:)
!      write(Fname,"(A,I1,A)") trim(Out_Dir)//"Dr_fcs_Logistic_based_on_dridx_for_"//trim(Dr_type)//"_len_",ilen,".nc"
!      Call Write_Dr_and_Verification (Fname,Dr_type,ilen)
      !---------------------------------------------------------------------
      CONTAINS
      !------------------------------------------------------------------------------
      Subroutine Write_Dr_and_Verification (Fname,Dr_type,ilen)
         Character (Len=*) :: Fname,Dr_type
         Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
         Integer :: ilead,ilen, it
         Real :: xlead(Nlead),xtime(Ntime)

         t_unit="months since 1983-01-01 "; var_unit="-";
         lon_nm="lon"; lat_nm="lat"; z_nm="lead"; t_nm="mon"
         Do ilead=1,Nlead; xlead(ilead)=ilead;  Enddo
         Do it=1,Ntime; xtime(it)=it-1; Enddo

         write(Var_ln,"(A,I2.2,A)") "Drought events of the ",ilen,"-month observed "//trim(Dr_type)//" drought"
         Call wr_ncdf_4d (Fname, Lon,Lat,xlead,xtime,Nlon,Nlat,Nlead,Ntime,Dr_obs,"obs", &
                          t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)
         write(Var_ln,"(A,I2.2,A)") "Drought events of the ",ilen,"-month forecast "//trim(Dr_type)//" drought"
         Call wr_ncdf_4d_apd (Fname,Dr_fcs,"fcs",Nlon,Nlat,Nlead,Ntime, &
                              Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss)

         Var_ln="Proportion Correction (PC): Time=1-6 corresponding to LT=1-6"
         Call wr_ncdf_3d_apd (Fname,PC,"pc",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss)
         Var_ln="Probability Of Detection (POD): Time=1-6 corresponding to LT=1-6"
         Call wr_ncdf_3d_apd (Fname,POD,"pod",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss)
         Var_ln="Threat Score (TS), or Critical Success Index (CSI): Time=1-6 corresponding to LT=1-6"
         Call wr_ncdf_3d_apd (Fname,TS,"ts",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss)
         Var_ln="Heidke Skill Score (HSS): Time=1-6 corresponding to LT=1-6"
         Call wr_ncdf_3d_apd (Fname,HSS,"hss",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss)
      End Subroutine Write_Dr_and_Verification
      !------------------------------------------------------------------------------
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"

      !------------------------------------------------------------------------------
      END PROGRAM t_logistic