!
!  CHUONG TRINH Du bao So luong XTND bang cac Mo hinh Thong ke
!  Cac Mo hinh thong ke duoc su dung:
!  1) ANN: 
!  2) HQTT: Hoi quy tuyen tinh nhieu bien (Ket qua giong voi HQTB)
!  3) LAD: Least Absolute Deviation
!  4) LMV: Minimax
!  Cac nhan to du bao: La cac EOF (3 EOF dau tien) cua 14 ClimIDX. So nhan to du
!  tuyen = 3 thang truoc thang du bao x 3 EOF = 9 nhan to
!  Cac nhan to du bao thuc te se duoc chon tu 9 nhan to du tuyen boi HQTB

   PROGRAM TC_Fcs
   Use Netcdf

!   Integer, Parameter :: Yr1=2012, Yr2=2017, Nmon=12, Nfile=Nmon, Nmodels=9
   Integer, Parameter :: Yr1=2003, Yr2=2009, Nmon=12, Nfile=Nmon, Nmodels=9
   Integer, Parameter :: NMax=Yr2-Yr1+1, MMax=10

   !Integer, Parameter :: NMax=1000, MMax=20
   Integer            :: Nyr1
   Real, Parameter    :: Miss=-99.

   REAL, ALLOCATABLE :: X(:,:), X1(:)              ! Input Data Arrays
   INTEGER, ALLOCATABLE :: IDX(:)                  ! Index of variables for stepwise regression

   Character (Len=70) :: FINP !, FOUT                ! Filename for input and output

   Real    :: Xorg(NMAX,MMax), Yorg(NMAX), Yhq(NMAX)
   ! 
   Character (Len=8) :: Ten(MMAX), Varname(MMax)        ! Variable's name
   Real :: Beta_tmp(MMax+1)                 ! Working array


   Character (Len=70) :: TIEUDE
   Character (Len=2)  :: ST1,ST2
   Character (Len=20) :: FMT

   !--------------------------------------------------------------------------------

   Integer :: ifile                    ! So file thuc te
   !Character (Len=100), Allocatable :: FnameInp(:), FnameOut(:), FnameSave(:)
   Integer :: Mvars, N, Mbien              ! So bien doc lap, Dung luong mau va Tat ca cac Bien (Mvars+1)

   !Character (Len=5), Allocatable :: Mdel_nm(:)
   Real, Allocatable :: Beta(:,:,:)   ! Cac he so HQ cua 3 mo hinh HQTT, LAD, MINIMAX,
   Real, Allocatable :: Yout(:,:,:)   ! Outputs of Yorg and Y-estimated by 4 models
   Real, Allocatable, Dimension(:,:) :: ME, MAE, RMSE, RME, Corre, Pvalue  ! Verification

   Integer :: Nyear, imon      ! Tong so nam co SL va chi so thang
   !--------------------------------------------------------------------------------
   ! ANN:
   character (len=250) :: input_file, output_file, save_file
   logical :: training_mode

   ! Use for calculating scaling factors
   Integer :: Nsample, Nepoc

   Real,allocatable, Dimension(:,:) :: Inputs, Outputs
   Real,allocatable, Dimension(:)   :: min_in, max_in, min_out, max_out
   Real :: mean_error
   !-------------------------------------------------------------------------------
    integer ::  Nin, Nhidden, Nout, iin, ihid, iout
    integer :: i,j,k

    Real,allocatable, Dimension(:,:) :: weight_ih, weight_ho, moment_ih, moment_ho
    Real,allocatable, Dimension(:)   :: fin, fhid, fout, ftarget, delta_o, delta_h
   ! Training
    Real :: learning_rate,momentum, error
    logical :: debug
    Character (Len=250) :: Y_out_Fname_NC, Y_out_Fname_TXT, Beta_Fname, Veri_Fname_NC, Veri_Fname_TXT
 
    !namelist /info/ Nhidden, learning_rate,momentum, training_mode,Nepoc
    !Namelist /NumFile/ Nfile                           ! Number of file to be calculated
    !Namelist /FileName/ FnameInp, FnameOut, FnameSave  ! Input and output filenames
    !Namelist /Number_Of_Model/ Nmodels
    !Namelist /Model_Names/ Mdel_nm
    !Namelist /Model_Training_outputs/ Y_out_Fname_NC, Y_out_Fname_TXT, Beta_Fname, Veri_Fname_NC, Veri_Fname_TXT 
   !--------------------------------------------------------------------------------
   !Integer :: Kmon, Lag, Kmon_0
   !Namelist /Lag_Kmon/ Lag, Kmon

   Real, Allocatable :: Mean_Obs(:), Std_Obs(:)
   !Integer, Parameter :: Yr1=1982, Yr2=2009, Nmon=12, Nfile=Nmon, Nmodels=9
   
   Character (Len=5), Parameter :: Mdel_nm(Nmodels+1)=(/ &
             "Obs","MLR","LAD","LMV","ANN","E12","E123","E124","E14","E1234"/)

   Integer :: Yr_Sta(Nmon), Nyrs(Nmon), NVars(Nmon)   ! Nam bat dau, Dung luong mau, So bien DL 
   Character (Len=250) :: Dir, Dir_fc
   Character (Len=250), Parameter :: Fnames(Nfile)=(/"tc3m_01.txt","tc3m_02.txt","tc3m_03.txt","tc3m_04.txt", &
                                                     "tc3m_05.txt","tc3m_06.txt","tc3m_07.txt","tc3m_08.txt", &
                                                     "tc3m_09.txt","tc3m_10.txt","tc3m_11.txt","tc3m_12.txt"/)
   !Character (Len=250), Parameter :: Fnames(Nfile)=(/"tc6m_01.txt","tc6m_02.txt","tc6m_03.txt","tc6m_04.txt", &
   !                                                  "tc6m_05.txt","tc6m_06.txt","tc6m_07.txt","tc6m_08.txt", &
   !                                                  "tc6m_09.txt","tc6m_10.txt","tc6m_11.txt","tc6m_12.txt"/)
   Character (Len=250) :: FnameOut(Nfile), FnameSave(Nfile), FnameInp(Nfile)
   !------------------- For ANN -------------------------------------------------------------
     Nhidden       =     2        ! Number of hidden layers
     learning_rate =   0.5
     momentum      =     0
     training_mode =.false.       ! For training
     Nepoc         =  8000        ! Iterative training times
   !--------------------------------------------------------------------------------
     Dir="/home/tanpv/Progs/TCs/FCS/Out/"
     Dir_fc="/home/tanpv/Progs/TCs/FCS/Out_Fcs/"
     Do i=1,Nfile 
        FnameInp(i)=trim(Dir_fc)//trim(Fnames(i)); FnameOut(i)=trim(Dir)//"Out_"//trim(Fnames(i)); FnameSave(i)=trim(Dir)//"Sav_"//trim(Fnames(i))
     Enddo
     Y_out_Fname_NC   = "/home/tanpv/Progs/TCs/FCS/Dat_out/Yout_fc.nc"
     Y_out_Fname_TXT  = "/home/tanpv/Progs/TCs/FCS/Dat_out/TC_out_fc.txt"
     Beta_Fname       = "/home/tanpv/Progs/TCs/FCS/Dat_out/Beta.nc"
     Veri_Fname_NC    = "/home/tanpv/Progs/TCs/FCS/Dat_out/Veri_Results_fc.nc"
     Veri_Fname_TXT   = "/home/tanpv/Progs/TCs/FCS/Dat_out/Veri_Results_fc.txt"
   !--------------------------------------------------------------------------------


   !Open (1, file="NameList_PCA.txt", Status="Old")
   !Read (1, Lag_Kmon)
   !Kmon_0 = Kmon + Lag + 1
   !Read (1, Number_Of_Model)
   !Allocate (Mdel_nm(Nmodels+1)); 
   !Read (1, Model_Names)
   !Read (1, Model_Training_outputs)
   !Close (1)
     

   !Open (1, file="NameList_TC_Training.txt")   ! Namelist file
   !Read (1,info)     ! For ANN
   !Read (1, NumFile); Allocate (FnameInp(Nfile), FnameOut(Nfile), FnameSave(Nfile)); Read (1, FileName)
   !Close (1)

     Allocate (Mean_Obs(Nfile), Std_Obs(Nfile))
     Allocate (Yout(Nfile,Nmax,0:Nmodels))
     Allocate (ME(Nfile,Nmodels),MAE(Nfile,Nmodels),RMSE(Nfile,Nmodels))
     Allocate (RME(Nfile,Nmodels),Corre(Nfile,Nmodels),Pvalue(Nfile,Nmodels))     

   Nyear=0; Yout(:,:,:) = Miss;
   Call Read_Beta (trim(Beta_Fname))

   Do ifile=1, Nfile

      Call Read_Data (ifile, FnameInp(ifile)) 
      imon = ifile; Mvars = NVars(imon)-1; Mbien=NVars(imon); N=Nyrs(imon)

      Yout(ifile,1:N,0) = Yorg(1:N) 
      Mean_Obs(ifile) = Ave  (Yorg(1:N), N, Miss)
      Std_Obs(ifile)  = Stdev(Yorg(1:N), N, Miss) 
 
      print*,"Calling Calc_HQTT..."
      Call Calc_HQTT (ifile)        ! HQTT
      Yout(ifile,1:N,1) = Yhq(1:N)
      print*,"Calling Calc_LAD..."

      Call Calc_LAD (ifile)        ! Least Absolute Deviation
      Yout(ifile,1:N,2) = Yhq(1:N)
      print*,"Calling Calc_LMV..."

      Call Calc_LMV (ifile)        ! Minimax
      Yout(ifile,1:N,3) = Yhq(1:N)
     !--------------------------ANN-------------------------
      debug = .false.
      input_file  = trim(FnameInp(ifile))
      output_file = trim(FnameOut(ifile))
      save_file   = trim(FnameSave(ifile))

      print*, "Calculating ",input_file
      Call controller_init

      if (training_mode) then       ! .T. = Training, .F. = forecast (Define from Namelist)
         Call controller_train_data
      endif

      write(*,*) "Writting output..."
      Call controller_output_data !(Yhq)
      Call controller_deinit
      Yout(ifile,1:N,4) = Yhq(1:N)
      !--------------------------ANN-------------------------

   Enddo

   ! Set all estimated negative TC-counts to Zero:
   Where (Yout/=Miss.AND.Yout<0) Yout=0

   Do ifile=1, Nfile
      ! Ensemble means: Yout(:,:,5) = (MLR+LAD+LMV)/3; Yout(:,:,6) = (MLR+LAD+LMV+ANN)/4
      ! "E12","E123","E124","E14","E1234"
      !   5     6      7      8      9
      ! Yr_Sta(Nmon), Nyrs(Nmon), NVars(Nmon)
      Nyear=Nyrs(ifile)
      Do iyr=1, Nyear
         if (Yout(ifile,iyr,1)/=Miss.AND.Yout(ifile,iyr,2)/=Miss) &
            Yout(ifile,iyr,5)=(Yout(ifile,iyr,1)+Yout(ifile,iyr,2))/2.
 
         if (Yout(ifile,iyr,1)/=Miss.AND.Yout(ifile,iyr,2)/=Miss.AND.Yout(ifile,iyr,3)/=Miss) &
            Yout(ifile,iyr,6)=(Yout(ifile,iyr,1)+Yout(ifile,iyr,2)+Yout(ifile,iyr,3))/3.

         if (Yout(ifile,iyr,1)/=Miss.AND.Yout(ifile,iyr,2)/=Miss.AND.Yout(ifile,iyr,4)/=Miss) &
            Yout(ifile,iyr,7)=(Yout(ifile,iyr,1)+Yout(ifile,iyr,2)+Yout(ifile,iyr,4))/3.

         if (Yout(ifile,iyr,1)/=Miss.AND.Yout(ifile,iyr,4)/=Miss) &
            Yout(ifile,iyr,8)=(Yout(ifile,iyr,1)+Yout(ifile,iyr,4))/2.
        
         if (Yout(ifile,iyr,1)/=Miss.AND.Yout(ifile,iyr,2)/=Miss.AND.Yout(ifile,iyr,3)/=Miss.AND.Yout(ifile,iyr,4)/=Miss) &
            Yout(ifile,iyr,9)=(Yout(ifile,iyr,1)+Yout(ifile,iyr,2)+Yout(ifile,iyr,3)+Yout(ifile,iyr,4))/4.
      Enddo
      ! Tinh sai so va HSTQ giua Yorg va Yhq
      Yorg(1:Nyear) = Yout(ifile,1:Nyear,0) 
      Do imod=1, Nmodels  
         print*,"Calling Verification of Model: ", imod
         Yhq(1:Nyear) = Yout(ifile,1:Nyear,imod) ! ; Where (Yhq/=Miss.AND.Yhq<0) Yhq=0
         Call Calc_Verification (Yorg, Yhq, Nyear, ME(ifile,imod), MAE(ifile,imod), RMSE(ifile,imod), RME(ifile,imod), Corre(ifile,imod), Pvalue(ifile,imod))
      Enddo
   Enddo

   print*,"Write_Yorg_Yhq_to_NC..."
   Call Write_Yorg_Yhq_to_NC (trim(Y_out_Fname_NC)) ! "Yout.nc")                      ! Gia tri TC observed va estimated

   print*,"Write Beta..."
   Call Write_Model_Coeffs_to_NC (trim(Beta_Fname)) !"Beta.nc")

   print*,"Writing to TXT file ..."
   Call Write_Yorg_Yhq_to_TXT (trim(Y_out_Fname_TXT)) !"TC_out.txt")

   print*,"Write_Verification_Results_to_NC"
   Call Write_Verification_Results_to_NC (trim(Veri_Fname_NC)) !"Veri_Results.nc")  ! Cac chi so danh gia
  
   Call Write_Verification_to_TXT (trim(Veri_Fname_TXT)) !"Veri_Results.txt")
   !--------------------------------------------------------------------------------
      
   CONTAINS
   !--------------------------------------------------------------------------------
   Subroutine Read_Beta (Fname)
     Character (Len=*) :: Fname
     Real, Allocatable :: Xlon(:), Xlat(:), Xtime(:)
     Integer :: Nfi, Ncoef

     Call rd_ncdf_3d (Fname,Xlon,Xlat,Xtime, Nfi, N_model, Ncoef, Beta, "Beta")
     print*, Nfi, N_model, Ncoef
     Deallocate (Xlon,Xlat,Xtime)
   End Subroutine Read_Beta
   !--------------------------------------------------------------------------------

Subroutine Write_Verification_to_TXT (Fname)
  Integer, Parameter :: Nmon=12
  Character (Len=*) :: Fname
  Character (Len=50) :: Region ! (/"VES_N15","VES_S15"/) 
  Integer :: imod, ifile, imon
  
  Open (3,file=trim(Fname), Status="Unknown")
  Write (3,*) "Verification results for the TC forecast using dependent datasets"
  Do ifile=1, Nfile
    ! if (ifile<=Nmon) then; Region="VES_N15"; else; Region="VES_S15"; endif
     Region="VES" !"trim(FnameInp(ifile))
     imon=ifile !mod(ifile,Nmon); if (imon==0) imon=Nmon
     Write (3,"(A,A,A,I2.2)") " Input file: ", trim(Region), " Mon: ",imon
     Write (3,"(2(A,F6.1))") " Mean TC-count: ", Mean_Obs(ifile), " Standard Deviation: ", Std_Obs(ifile)
     Write (3,"(A8,20A8)") "ID/Model", (Mdel_nm(imod),imod=2,Nmodels+1)
     Write (3,"(A8,20F8.1)") " ME: ", (ME(ifile,imod), imod=1, Nmodels) 
     Write (3,"(A8,20F8.1)") " MAE: ", (MAE(ifile,imod), imod=1, Nmodels)
     Write (3,"(A8,20F8.1)") " RMSE: ", (RMSE(ifile,imod), imod=1, Nmodels)
     Write (3,"(A8,20F8.1)") " RMAE: ", (RME(ifile,imod), imod=1, Nmodels)
     Write (3,"(A8,20F8.3)") " CORR: ", (Corre(ifile,imod), imod=1, Nmodels)
     Write (3,"(A8,20F8.3)") " Pval: ", (Pvalue(ifile,imod), imod=1, Nmodels)
     Write (3,"(20A8)") ("---------", i=1,6)
  Enddo
End Subroutine Write_Verification_to_TXT
!--------------------------------------------------------------------------------

Subroutine Write_Model_Coeffs_to_NC (Fname)
  Character (Len=*) :: Fname
  Character (Len=255) :: var_ln, t_unit
  Real, Allocatable :: Xlon(:), Xlat(:), Xtime(:), X(:,:,:)
  Integer, Parameter :: NModel=3 ! MLR, LAD, LMV
  Integer :: Mvars

  Mvars=Mbien-1
  Allocate (Xlon(Nfile), Xlat(NModel), Xtime(Mvars+1), X(Nfile,Nmodel,Mvars+1))
  Do ilon=1, Nfile; Xlon(ilon)=ilon; Enddo
  Do ilat=1, Nmodel; Xlat(ilat)=ilat; Enddo
  Do it=1, Mvars+1; Xtime(it)=it-1; Enddo
  Do ifile=1, Nfile
     Do imod=1, Nmodel
        X(ifile,imod,1:Mvars+1) = Beta(ifile,imod,0:Mvars)
     Enddo
  Enddo
  !print*,X
  print*,"Nfile, Nmodel, Mvars=", Nfile, Nmodel, Mvars+1
  t_unit="years since 1961 "
  var_ln="Cac he so cua cac PTHQ b0, b1,... cua cac mo hinh MLR, LAD, LMV"
  Call wr_ncdf_3d (trim(Fname),XLon,XLat,Xtime, Nfile,Nmodel,Mvars+1, X, "Beta", & 
                   trim(t_unit), trim(var_ln),"-","files","Model","heso",Miss)

End Subroutine Write_Model_Coeffs_to_NC
!--------------------------------------------------------------------------------
Subroutine Write_Verification_Results_to_NC (Fname)
  Character (Len=255) :: var_ln
  Character (Len=*) :: Fname
  Real, Allocatable :: Xlon(:), Xlat(:), X(:,:)
  Integer :: ilon, ilat, it

  Allocate (Xlon(Nfile),Xlat(Nmodels), X(Nfile, Nmodels))
  Do ilon=1, Nfile; Xlon(ilon) = ilon; Enddo      ! So file
  Do ilat=1, Nmodels; Xlat(ilat) = ilat; Enddo    ! 4 Models

  var_ln="Mean Error (ME) forecast of TC-count by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = ME(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d     (trim(Fname), Xlon, Xlat, Nfile, Nmodels, X, "ME", "Mean Error", "NTC", "Case", "Model", Miss)

  var_ln="Mean Absolute Error (MAE) forecast of TC-count by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = MAE(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d_apd (trim(Fname), X, "MAE", Nfile, Nmodels, "Mean Absolute Error", "NTC", "Case", "Model", Miss)

  var_ln="Root Mean Square Error (RMSE) forecast of TC-count by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = RMSE(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d_apd (trim(Fname), X, "RMSE", Nfile, Nmodels, "Root Mean Square Error", "NTC", "Case", "Model", Miss)

  var_ln="Mean Relative Absolute Error (RME) forecast of TC-count by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = RME(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d_apd (trim(Fname), X, "RME", Nfile, Nmodels, "Mean Relative Absolute Error", "%", "Case", "Model", Miss)

  var_ln="Multi Correlation Coefficients (Corr) of TC-count forecast by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = Corre(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d_apd (trim(Fname), X, "Corr", Nfile, Nmodels, "Multi Correlation Coefficients", "-", "Case", "Model", Miss)

  var_ln="P-value of Multi Correlation Coefficients (Corr) of TC-count forecast by different models for Jan to Dec and for North and South VES (Case)"
  X(1:Nfile,1:Nmodels) = Pvalue(1:Nfile,1:Nmodels)
  Call wr_ncdf_2d_apd (trim(Fname), X, "Prob", Nfile, Nmodels, "P-value of Multi Correlation Coefficients", "-", "Case", "Model", Miss)
  Deallocate (Xlon, Xlat, X)
End Subroutine Write_Verification_Results_to_NC
!--------------------------------------------------------------------------------
Subroutine Write_Yorg_Yhq_to_TXT (Fname)
  Character (Len=*) :: Fname
  Real, Allocatable :: X(:,:,:) ! (Nfile, Nmodels+1, N)
  Integer :: ifile, iyr, imon, imod, it1, it2, Nt
  ! Yr_Sta(Nmon), Nyrs(Nmon)
  
  Nyear=Nmax
  Open (3,file=trim(Fname),status="Unknown")
  Allocate (X(1:Nfile,1:Nmodels+1,1:Nyear)); X(:,:,:) = Miss

  Do ifile=1, Nfile

     write(3,*) "Input file: ",trim(FnameInp(ifile))
     write(3,"(A5,10A7)") "Yrs", (Mdel_nm(imod),imod=1,Nmodels+1)
     imon=ifile
     it1=Yr_Sta(imon)-Yr1+1; it2=it1+Nyrs(imon)-1; Nt=Nyrs(imon)
     Do imod=1, Nmodels+1
        !X(ifile,imod,1:Nyear) = Yout(ifile,1:Nyear,imod-1)
        X(ifile,imod,it1:it2) = Yout(ifile,1:Nt,imod-1)
     Enddo

     Do iyr=1, Nyear
        write(3,"(I5,10F7.1)") iyr+Yr1-1, (X(ifile,imod,iyr), imod=1,Nmodels+1)
     Enddo
  Enddo
  Deallocate (X); Close (3)
End Subroutine Write_Yorg_Yhq_to_TXT
!--------------------------------------------------------------------------------
Subroutine Write_Yorg_Yhq_to_NC (Fname)
  Character (Len=255) :: t_unit, var_ln
  Character (Len=*) :: Fname
  Real, Allocatable :: X(:,:,:) ! (Nfile, Nmodels+1, N)
  Real, Allocatable :: Xlon(:), Xlat(:), Xtime(:)
  Integer :: ilon, ilat, it, imon, it1, it2, Nt

  ! Yr_Sta(Nmon), Nyrs(Nmon)
  Nyear=Nmax
  Allocate (X(1:Nfile,1:Nmodels+1,1:Nyear),Xlon(Nfile),Xlat(Nmodels+1),Xtime(Nyear)); X(:,:,:) = Miss
  Do ilon=1, Nfile; Xlon(ilon)=ilon; Enddo     ! So truong hop du bao: 12 thang x 2 vung = 24
  Do ilat=1, Nmodels+1; Xlat(ilat)=ilat; Enddo ! Yorg + 4 models
  Do it=1, Nyear; Xtime(it)=Yr1+it-1; Enddo              ! So nam du bao = 1962..2010
  write (t_unit,"(A,I4.4,A)") "years since ",Yr1,"-01-01 "
  var_ln="Observed TC-count and its Estimation by different models (MLR, LAD, LMV, ANN)"
  X=Miss
  Do ifile=1, Nfile
     imon=ifile !mod(ifile,12); if (imon==0) imon=12 
     it1=Yr_Sta(imon)-Yr1+1; it2=it1+Nyrs(imon)-1; Nt=Nyrs(imon) 
     Do imod=1, Nmodels+1
        !X(ifile,imod,1:Nyear) = Yout(ifile,1:Nyear,imod-1)
        X(ifile,imod,it1:it2) = Yout(ifile,1:Nt,imod-1)
     Enddo
  Enddo
  Call wr_ncdf_3d(trim(Fname),XLon,XLat,Xtime, Nfile,Nmodels+1,Nyear, X, "tc", trim(t_unit), trim(var_ln),"1","mon","model","year",Miss) 
  ! wr_ncdf_3d (filenm,lon,lat,t,nlon,nlat,nt,var,var_nm,t_unit, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag)
  Deallocate (Xlon, Xlat, Xtime, X)
End Subroutine Write_Yorg_Yhq_to_NC 
!--------------------------------------------------------------------------------

Subroutine Calc_HQTT (ifile)
  Beta_tmp(1:Mvars+1) = Beta(ifile,1,1:Mvars+1)
  print*,Beta_tmp(1:Mvars+1)
  Do i=1, N
     Yhq(i) = 0
     Do j=2, Mvars+1
        Yhq(i) = Yhq(i)+Beta_tmp(j)*Xorg(i,j-1) ! Beta_tmp(2..Mvars+1)=A1,...Am
     Enddo
     Yhq(i) = Yhq(i)+Beta_tmp(1)                ! Beta_tmp(1)=A0
  Enddo
End Subroutine Calc_HQTT
!--------------------------------------------------------------------------------
Subroutine Calc_LAD (ifile)
  Integer :: i, j, IOPT

  Beta_tmp(1:Mvars+1) = Beta(ifile,2,1:Mvars+1)
  Do i=1, N
     Yhq(i) = 0
     Do j=2, Mvars+1
        Yhq(i) = Yhq(i)+Beta_tmp(j)*Xorg(i,j-1) ! Beta_tmp(2..Mvars+1)=A1,...Am
     Enddo
     Yhq(i) = Yhq(i)+Beta_tmp(1)                ! Beta_tmp(1)=A0
  Enddo
End Subroutine Calc_LAD
!--------------------------------------------------------------------------------
Subroutine Calc_LMV (ifile)
  Integer :: i, j, IOPT

  Beta_tmp(1:Mvars+1) = Beta(ifile,3,1:Mvars+1)
  Do i=1, N
     Yhq(i) = 0
     Do j=2, Mvars+1
        Yhq(i) = Yhq(i)+Beta_tmp(j)*Xorg(i,j-1) ! Beta_tmp(2..Mvars+1)=A1,...Am
     Enddo
     Yhq(i) = Yhq(i)+Beta_tmp(1)                ! Beta_tmp(1)=A0
  Enddo
End Subroutine Calc_LMV
!--------------------------------------------------------------------------------

Subroutine Calc_Verification (Yorg,Yhq,N, ME, MAE, RMSE, RME, Corr, Pvalue)
! Input: Yorg(1:N) is an array stored the origin values of Y
!        Xorg(1:N,1:M) is an array stored the origin values of X
!        Beta(0:M) is an array stored the regression coeficients B0, B1,..., BM
! Output: ME is Mean Error
!         MAE is Mean Absolute Error
!         RMSE is Root Mean Square Error
!         RME is Relative Mean Error
!         Corr, Pvalue are Correlation Coefficient and P-value (t-test)
Integer :: N, i, j, k, Nc, Nc0
Real    :: Yorg(N), Yhq(N), Error
Real    :: ME, MAE, RMSE, RME, Corr, Pvalue

  ME=0; MAE=0; RMSE=0; RME=0; Nc=0; Nc0=0 
  Do i=1, N
     if (Yhq(i)/=Miss.AND.Yorg(i)/=Miss) then
        Nc=Nc+1
        Error= Yhq(i)-Yorg(i)
        ME   = ME   + Error
        MAE  = MAE  + ABS(Error)
        RMSE = RMSE + Error*Error
        if (Yorg(i)/=0) then; Nc0=Nc0+1; RME  = RME  + ABS(Error/Yorg(i)); endif
     endif
  Enddo
  ME = ME/Real(Nc); MAE = MAE/Real(Nc); RME = RME/Real(Nc0)*100.; RMSE = Sqrt(RMSE/Real(Nc))
  Call Correlation (Yorg, Yhq, N, Corr, Pvalue, Miss)
End Subroutine Calc_Verification
!--------------------------------------------------------------------------------
Subroutine Read_Data (ifile, Finp) ! Calc_Stepwise_Regression (ifile, Finp, Fout)
   Character (Len=*) :: Finp !, Fout
   Integer :: ifile, iv_y, Nc, Nyr1   ! Chi so bien phu thuoc
   Real, Allocatable :: Xtmp(:)
   Real              :: Tmp

      OPEN(1,FILE=trim(FINP),Status="Old")

      PRINT*," READ IN DATA FROM FILE "//trim(FINP) !,"  Output file:",trim(FOUT)

      READ (1,"(A)") TIEUDE
      Read(1,*) Nyr1, N, Mbien
      ALLOCATE (IDX(Mbien), Xtmp(Mbien-1))

      Read(1,*) (Idx(i),i=1,Mbien)
      Read (1,"(20A10)") (Ten(i),i=1,Mbien)
      Nc=0
   L1:   DO I=1,N
         READ(1,*) Tmp, (Xtmp(j), j=1,Mbien-1)
         if (Tmp==Miss) Cycle L1
         Do j=1, Mbien-1
            if (Xtmp(j)==Miss) Cycle L1
         Enddo
         Nc=Nc+1;
         Yorg(Nc)=Tmp; Xorg(Nc,1:Mbien-1)=Xtmp(1:Mbien-1)
!         READ(1,*) Yorg(i), (Xorg(i,j), j=1,Mbien-1) ! (X(I,J),J=1,Mbien)
      ENDDO L1
      N = Nc
      CLOSE(1)
      Deallocate (IDX,Xtmp)
      !
      Yr_Sta(ifile)=Nyr1; Nyrs(ifile)=N; NVars(ifile)=Mbien

End Subroutine Read_Data 
!--------------------------------------------------------------------------------
     
Include "/home/tanpv/Mylibs/Sub_HQTT.f90"
Include "/home/tanpv/Mylibs/Sub_LAD.f90"
Include "/home/tanpv/Mylibs/Sub_LMV.f90"
Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90"
Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" 
Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
Include "/home/tanpv/Mylibs/Sub_neural.f90"
Include "/home/tanpv/Mylibs/Stdev.f90"
Include "/home/tanpv/Mylibs/Ave.f90"
!--------------------------------------------------------------------------------
END
Include "/home/tanpv/Mylibs/Student_Prob.f90"
!-------------------------------------------------------------------------------------------