program main

!*****************************************************************************80
!! MAIN is the main program for KMEANS_TEST.
!  Discussion:
!    KMEANS_TEST tests the KMEANS library.
!  Licensing:
!    This code is distributed under the GNU LGPL license.
!  Modified:
!    06 October 2009
!  Author:
!    John Burkardt
!
  Use NetCDF
  implicit none
  Real, Parameter :: Rmiss=-99.0
  Integer, Parameter :: Imiss=-99
  Real, Allocatable :: StnID(:), X(:,:), Lon(:), Lat(:), Reg(:), xtime(:)
  Integer :: Nst,Ntime,i,j, irow, icol, N_Clus

  Real (kind=8), Allocatable :: point(:,:)   ! (M,N)
  Real, Allocatable :: C(:,:), XX(:,:)
  Integer,Allocatable :: Z(:)

  Integer (kind=4) :: M, N, cluster_num  ! Dims, Points, Num of clusters
  Logical :: status

  Integer, Allocatable :: Clus(:,:)      ! (N,N_Clus)
  Real,    Allocatable :: Clus_var(:,:)  ! (N_Clus,N_Clus)
  Real,    Allocatable :: Total_var(:)   ! (N_Clus)
  Real,    Allocatable :: Total_Ene(:)   ! (N_Clus)
  Integer, Allocatable :: Clus_popu(:,:) ! (N_Clus,N_Clus)

  integer ( kind = 4 ), allocatable :: cluster_o(:)
  real ( kind = 8 ), allocatable :: cluster_center_o(:,:)
  real ( kind = 8 ), allocatable :: cluster_energy_o(:)
  integer ( kind = 4 ), allocatable :: cluster_population_o(:)
  real ( kind = 8 ), allocatable :: cluster_variance_o(:)

  !allocate ( cluster(point_num) )
  !allocate ( cluster_population(cluster_num) )
  !allocate ( cluster_variance(cluster_num) )

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'KMEANS_TEST'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Test the KMEANS library.'

  ! TEST08 applies KMEANS_01 + KMEANS_03 to points_100.txt:
  Call Read_data          ! Return: X(Nst,Ntime), Lon(Nst),Lat(Nst),Reg(Nst)

  M=Ntime; N=Nst

  Allocate (point(M,N), XX(M,N))
  Do irow=1,N
     point(1:M,irow) = X(irow,1:M)
     XX(1:M,irow) = X(irow,1:M)
  Enddo

  !cluster_num=1
  !Call test03 (point,M,N, cluster_num, Reg)
  !pause

  !cluster_num=10
  N_Clus=10
  Allocate (Clus(Nst,N_Clus),Clus_var(N_Clus,N_Clus),Clus_popu(N_Clus,N_Clus))
  Allocate (Total_Ene(N_Clus),Total_var(N_Clus))
  Clus=Imiss; Clus_popu=Imiss; Clus_var=Rmiss; Total_var=Rmiss

  Do j=1,N_Clus
     cluster_num=j
     allocate ( cluster_o(N) )
     allocate ( cluster_center_o(M,N) )
     allocate ( cluster_energy_o(cluster_num) )
     allocate ( cluster_population_o(cluster_num) )
     allocate ( cluster_variance_o(cluster_num) )

     !call test04 (point,M,N,cluster_num)
     !Call test03 (point,M,N, cluster_num, Reg)
     Call test08 (point,M,N, cluster_num, Reg)
     print*,"Case ",j

     Clus(1:N,j)=cluster_o(1:N)
     Clus_var(1:j,j)=cluster_variance_o(1:cluster_num)
     Clus_popu(1:j,j)=cluster_population_o(1:cluster_num)
     Total_var(j)=sum(cluster_variance_o(1:cluster_num))/cluster_num
     Total_Ene(j)=sum(cluster_energy_o(1:cluster_num))

     Deallocate (cluster_o,cluster_center_o,cluster_energy_o,cluster_population_o,cluster_variance_o)
  Enddo

  pause
 ! Clusters:
  print"(2A4,20(3x,A1,I2.2))","STT","Reg",(("C",i),i=2,N_Clus)
  Do i=1,Nst
     print"(2I4,20I6)",i,Int(Reg(i)),(Clus(i,2:N_Clus))
  Enddo

 !Populations:
  print"(A3,20(3x,A1,I2.2))","Gr",(("C",i),i=2,N_Clus)
  Do i=1,N_Clus
     print"(I2,20I6)", i, (Clus_popu(i,2:N_Clus))
  Enddo

 ! Variances:
  print"(A3,20(6x,A1,I2.2))","Gr",(("C",i),i=2,N_Clus)
  Do i=1,N_Clus
     print"(I2,20F9.2)", i, (Clus_var(i,2:N_Clus))
  Enddo

 ! Total variances:
  print"(A3,2A12)","Gr","Tot_Var","Tot_Ene"
  Do i=2,N_Clus
     print"(I2,20F12.2)", i, Total_var(i), Total_Ene(i)
  Enddo

!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'KMEANS_TEST'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  ! Another method:
  Allocate (C(M,N_Clus),Z(N))
  Call KMPP (XX, M, N, N_Clus, C, Z)
  Do i=1,N
     print*,i,Int(Reg(i)),Z(i)
  Enddo
  pause

  Deallocate (StnID,X,Lon,Lat,Reg,xtime,point,XX, C, Z)

  stop 0
  !-----------------------------------------------------------------------------
  CONTAINS
      !--------------------------------------------------------------------
      Subroutine Read_data
        Character (Len=200) :: Fname

        Fname="/home/tanpv/SEACLID_New/RES/PDSI_Res/"// &
              "PDSI_Stn_No_Rain_Center.nc"
        Call rd_ncdf_2d (Fname,StnID,Xtime,Nst,Ntime,X,"ori_pdsi")
        !Call rd_ncdf_2d (Fname,StnID,Xtime,Nst,Ntime,X,"Pre")
        Deallocate (StnID)
        Call rd_ncdf_1d (Fname, StnID, Nst, Lon, "lon")
        Deallocate (StnID)
        Call rd_ncdf_1d (Fname, StnID, Nst, Lat, "lat")
        Deallocate (StnID)
        Call rd_ncdf_1d (Fname, StnID, Nst, Reg, "reg")
        print*, Reg
      End Subroutine Read_data

  Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
!!!end
!-----------------------------------------------------------------------------
 Include "kmean_lib.f90"
 Include "Sub_kmean.f90"
!-----------------------------------------------------------------------------
subroutine test04 (point,dim_num,point_num, cluster_num) !, &
           !cluster,cluster_population,cluster_variance,cluster_energy)

!  Clus(1:point_num,cluster_num)=cluster(1:point_num)
!  Clus_var(1:cluster_num,cluster_num)=cluster_variance(1:cluster_num)
!  Clus_popu(1:cluster_num,cluster_num)=cluster_population(1:cluster_num)
!  Total_var(cluster_num)=sum(cluster_variance(1:cluster_num))/cluster_num
!  Total_Ene(cluster_num)=sum(cluster_energy(1:cluster_num))

!*****************************************************************************80
!
!! TEST04 tries out the KMEANS_02 routine.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    06 October 2009
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ) cluster_num
  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  character ( len = 255 ) :: center_filename = 'test04_centers.txt'
  integer ( kind = 4 ), allocatable :: cluster(:)
  real ( kind = 8 ), allocatable :: cluster_center(:,:)
  real ( kind = 8 ), allocatable :: cluster_energy(:)
  character ( len = 255 ) :: cluster_filename = 'test04_clusters.txt'
  integer ( kind = 4 ), allocatable :: cluster_population(:)
  real ( kind = 8 ), allocatable :: cluster_variance(:)
  integer ( kind = 4 ) it_max
  integer ( kind = 4 ) it_num
  !real ( kind = 8 ), allocatable, dimension ( :, : ) :: point

  character ( len = 255 ) :: point_filename = 'points_100.txt'
  integer ( kind = 4 ) seed

  Integer :: i,j
  real ( kind = 8 ) :: point(dim_num,point_num)
  !Real ::  Reg(point_num)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  Test the KMEANS_02 algorithm.'
  write ( *, '(a)' ) '  (Applied Statistics Algorithm #136)'
!
!  Read the data points.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Data points will be read from "' // trim ( point_filename ) // '".'

  !call r8mat_header_read ( point_filename, dim_num, point_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Point spatial dimension = ', dim_num
  write ( *, '(a,i8)' ) '  Number of points = ', point_num

  !allocate ( point(dim_num,point_num) )

  !call r8mat_data_read ( point_filename, dim_num, point_num, point )
!
!  Clustering parameters.
!
  !cluster_num = 1
  it_max = 50
  seed = 123456789

  allocate ( cluster(point_num) )
  allocate ( cluster_center(dim_num,cluster_num) )
  allocate ( cluster_energy(cluster_num) )
  allocate ( cluster_population(cluster_num) )
  allocate ( cluster_variance(cluster_num) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of iterations allowed is ', it_max
!
!  Initialize the centers.
!
  print*,"initializing..."
  call cluster_initialize_1 ( dim_num, point_num, cluster_num, point, &
    cluster_center )

  cluster(1:cluster_num) = 0

  print*,"kmeans_02..."
  call kmeans_02 ( dim_num, point_num, cluster_num, it_max, it_num, point, &
    cluster, cluster_center, cluster_population, cluster_energy )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of iterations taken is ', it_num

  print*,"variance_compute..."
  call cluster_variance_compute ( dim_num, point_num, cluster_num, point, &
    cluster, cluster_center, cluster_variance )

  !call cluster_print_summary ( point_num, cluster_num, &
  !  cluster_population, cluster_energy, cluster_variance )

  !call r8mat_write ( center_filename, dim_num, cluster_num, &
  !  cluster_center )

  print*,">>>>>>"
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) &
    '  Cluster centers written to "' // trim ( center_filename ) // '".'

  !call i4mat_write ( cluster_filename, 1, point_num, cluster )

  !Clus(1:point_num,cluster_num)=cluster(1:point_num)
  !Clus_var(1:cluster_num,cluster_num)=cluster_variance(1:cluster_num)
  !Clus_popu(1:cluster_num,cluster_num)=cluster_population(1:cluster_num)
  !Total_var(cluster_num)=sum(cluster_variance(1:cluster_num))/cluster_num
  !Total_Ene(cluster_num)=sum(cluster_energy(1:cluster_num))

  !cluster_o(1:point_num)=cluster(1:point_num)
  !cluster_variance_o(1:cluster_num)=cluster_variance(1:cluster_num)
  !cluster_population_o(1:cluster_num)=cluster_population(1:cluster_num)
  !cluster_energy_o(1:cluster_num)=cluster_energy(1:cluster_num)

  cluster_o            = cluster
  cluster_variance_o   = cluster_variance
  cluster_population_o = cluster_population
  cluster_energy_o     = cluster_energy


  !Do i=1,point_num
  !   print*,i,Int(Reg(i)),cluster(i)
  !Enddo

  write ( *, '(a)' ) &
    '  Cluster assignments written to "' // trim ( cluster_filename ) // '".'

  deallocate ( cluster )
  deallocate ( cluster_center )
  deallocate ( cluster_energy )
  deallocate ( cluster_population )
  deallocate ( cluster_variance )
  !deallocate ( point )

  return
end
!*****************************************************************************80

subroutine test03 (point,dim_num,point_num, cluster_num, Reg)

!*****************************************************************************80
!
!! TEST03 tries out the KMEANS_01 routine.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    06 October 2009
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ) cluster_num
  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  character ( len = 255 ) :: center_filename = 'test03_centers.txt'
  integer ( kind = 4 ), allocatable :: cluster(:)
  real ( kind = 8 ), allocatable :: cluster_center(:,:)
  real ( kind = 8 ), allocatable :: cluster_energy(:)
  character ( len = 255 ) :: cluster_filename = 'test03_clusters.txt'
  integer ( kind = 4 ), allocatable :: cluster_population(:)
  real ( kind = 8 ), allocatable :: cluster_variance(:)
  integer ( kind = 4 ) it_max
  integer ( kind = 4 ) it_num
  !real ( kind = 8 ), allocatable, dimension ( :, : ) :: point
  character ( len = 255 ) :: point_filename = 'points_100.txt'
  integer ( kind = 4 ) seed

  Integer :: i,j
  real ( kind = 8 ) :: point(dim_num,point_num)
  Real ::  Reg(point_num)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  Test the KMEANS_01 algorithm.'
  write ( *, '(a)' ) '  (Applied Statistics Algorithm #58)'
!
!  Read the data points.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Data points will be read from "' // trim ( point_filename ) // '".'

  !call r8mat_header_read ( point_filename, dim_num, point_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Point spatial dimension = ', dim_num
  write ( *, '(a,i8)' ) '  Number of points = ', point_num

  !allocate ( point(dim_num,point_num) )

  !call r8mat_data_read ( point_filename, dim_num, point_num, point )
!
!  Clustering parameters.
!
  !cluster_num = 5
  it_max = 50
  seed = 123456789

  allocate ( cluster(point_num) )
  allocate ( cluster_center(dim_num,cluster_num) )
  allocate ( cluster_energy(cluster_num) )
  allocate ( cluster_population(cluster_num) )
  allocate ( cluster_variance(cluster_num) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of iterations allowed is ', it_max
!
!  Initialize the centers.
!
  call cluster_initialize_5 ( dim_num, point_num, cluster_num, point, &
    seed, cluster_center )

  cluster(1:cluster_num) = 0

  call kmeans_01 ( dim_num, point_num, cluster_num, it_max, it_num, point, &
    cluster, cluster_center, cluster_population, cluster_energy )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of KMEANS_01 iterations taken is ', it_num

  call cluster_variance_compute ( dim_num, point_num, cluster_num, point, &
    cluster, cluster_center, cluster_variance )

  call cluster_print_summary ( point_num, cluster_num, &
    cluster_population, cluster_energy, cluster_variance )

  call r8mat_write ( center_filename, dim_num, cluster_num, &
    cluster_center )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) &
    '  Cluster centers written to "' // trim ( center_filename ) // '".'

  !call i4mat_write ( cluster_filename, 1, point_num, cluster )

  Do i=1,point_num
     print*,i,Int(Reg(i)),cluster(i)
  Enddo

  cluster_o            = cluster
  cluster_variance_o   = cluster_variance
  cluster_population_o = cluster_population
  cluster_energy_o     = cluster_energy

  write ( *, '(a)' ) &
    '  Cluster assignments written to "' // trim ( cluster_filename ) // '".'

  deallocate ( cluster )
  deallocate ( cluster_center )
  deallocate ( cluster_energy )
  deallocate ( cluster_population )
  deallocate ( cluster_variance )
  !deallocate ( point )

  return
end

!*****************************************************************************80

subroutine test08 (point,dim_num,point_num, cluster_num, Reg)

!*****************************************************************************80
!
!! TEST08 tries out the HMEANS_01 + KMEANS_03 routine.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    06 October 2009
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer ( kind = 4 ) cluster_num
  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  character ( len = 255 ) :: center_filename = 'test08_centers.txt'
  integer ( kind = 4 ), allocatable :: cluster(:)
  real ( kind = 8 ), allocatable :: cluster_center(:,:)
  real ( kind = 8 ), allocatable :: cluster_energy(:)
  character ( len = 255 ) :: cluster_filename = 'test08_clusters.txt'
  integer ( kind = 4 ), allocatable :: cluster_population(:)
  real ( kind = 8 ), allocatable :: cluster_variance(:)
  integer ( kind = 4 ) it_max
  integer ( kind = 4 ) it_max_h
  integer ( kind = 4 ) it_max_k
  integer ( kind = 4 ) it_num
  Integer :: i,j

  real ( kind = 8 ) :: point(dim_num,point_num)
  Real ::  Reg(point_num)

  !real ( kind = 8 ), allocatable, dimension ( :, : ) :: point

  character ( len = 255 ) :: point_filename = 'points_100.txt'
  integer ( kind = 4 ) seed

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  Test the HMEANS_01 + KMEANS_03 algorithm.'
!
!  Read the data points.
!
  
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Data points will be read from "' // trim ( point_filename ) // '".'

  !call r8mat_header_read ( point_filename, dim_num, point_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Point spatial dimension = ', dim_num
  write ( *, '(a,i8)' ) '  Number of points = ', point_num

! dim_num   = Number of columns
! point_num = Number of points

  !allocate ( point(dim_num,point_num) )

  !call r8mat_data_read ( point_filename, dim_num, point_num, point )
!
!  Clustering parameters.
!
  !cluster_num = 3
  it_max = 50
  seed = 123456789
  it_max_h = 3
  it_max_k = it_max

  allocate ( cluster(point_num) )
  allocate ( cluster_center(dim_num,cluster_num) )
  allocate ( cluster_energy(cluster_num) )
  allocate ( cluster_population(cluster_num) )
  allocate ( cluster_variance(cluster_num) )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Initialize by using a few steps of HMEANS_02:'
  write ( *, '(a,i8)' ) '  Number of HMEANS_01 iterations allowed is ', it_max_h
  write ( *, '(a,i8)' ) '  Number of KMEANS_03 iterations allowed is ', it_max_k
!
!  Initialize the centers.
!
  call cluster_initialize_5 ( dim_num, point_num, cluster_num, point, &
    seed, cluster_center )
!
!  Initialize the clusters.
!
  cluster(1:cluster_num) = 0

  call hmeans_01 ( dim_num, point_num, cluster_num, it_max_h, it_num, point, &
    cluster, cluster_center, cluster_population, cluster_energy )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of HMEANS_01 iterations taken is ', it_num

  call cluster_variance_compute ( dim_num, point_num, cluster_num, point, &
    cluster, cluster_center, cluster_variance )

  call cluster_print_summary ( point_num, cluster_num, &
    cluster_population, cluster_energy, cluster_variance )

  call kmeans_03 ( dim_num, point_num, cluster_num, it_max_k, it_num, point, &
    cluster, cluster_center, cluster_population, cluster_energy )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of KMEANS_03 iterations taken is ', it_num

  call cluster_variance_compute ( dim_num, point_num, cluster_num, point, &
    cluster, cluster_center, cluster_variance )

  call cluster_print_summary ( point_num, cluster_num, &
    cluster_population, cluster_energy, cluster_variance )

  !call r8mat_write ( center_filename, dim_num, cluster_num, &
  !  cluster_center )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) &
    '  Cluster centers written to "' // trim ( center_filename ) // '".'

  !call i4mat_write ( cluster_filename, 1, point_num, cluster )

  cluster_o            = cluster
  cluster_variance_o   = cluster_variance
  cluster_population_o = cluster_population
  cluster_energy_o     = cluster_energy

  Do i=1,point_num
     print*,i,Int(Reg(i)),cluster(i)
  Enddo

  write ( *, '(a)' ) &
    '  Cluster assignments written to "' // trim ( cluster_filename ) // '".'

  deallocate ( cluster )
  deallocate ( cluster_center )
  deallocate ( cluster_energy )
  deallocate ( cluster_population )
  deallocate ( cluster_variance )
  !deallocate ( point )

  return
end
!------------------ End Prog. ------------------
END