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, Nmon=12
  Real, Allocatable :: X(:,:), Lon(:), Lat(:), Reg(:,:), Reg_pca(:,:), xtime(:)
  Real, Allocatable :: Mask(:,:), xlon(:), xlat(:)
  Integer :: Landpoint, Nst, ist, Nlat, Nlon, ilat, ilon, Ntime,i,j, irow, icol, N_Clus

  Real (kind=8), Allocatable :: point(:,:),pca(:,:)   ! (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(:)

  !--------------------------------------------------------------------------
  Integer, Parameter :: NClu_test=1 !4
  !Integer, Parameter :: Clus_test(NClu_test)=(/5,7,9,11/)
  !Integer, Parameter :: Clus_test(NClu_test)=(/5,9,11,15/)
  Integer, Parameter :: Clus_test(NClu_test)=(/5/)
  Integer :: iclu
  Character (Len=15) :: Dat_src="APH"
!  Character (Len=15) :: Dat_src="Ens"

  Character (Len=15) :: varn="spei_12"
!  Character (Len=15) :: varn="spi_12"

  Character (Len=255) :: Fname
  Logical :: PCA_Dat
  !--------------------------------------------------------------------------

    !Fname="/home/tanpv/Progs/Dr_SEA/Results/BL_"//trim(Dat_src)//".nc"

    !Fname="/home/tanpv/Progs/Dr_SEA/Res_APH/BL_"//trim(Dat_src)//".nc"
    !PCA_Dat=.False.

    Fname="/home/tanpv/Progs/Dr_SEA/EOF_New/pca_"//trim(Dat_src)//".nc"
    PCA_Dat=.True.; 
    if (PCA_Dat) varn="spe12"


    ! .True. if read in PCA, .False. if read in drought index
    Call Read_data(Fname,varn,PCA_Dat)          ! Return: X(Nst,Ntime), Lon(Nst),Lat(Nst), Mask(Nlon,Nlat)

    ! Number of clusters:
    N_Clus=13

    ! Calculate using drought index:
     M=Ntime; N=Landpoint !Nst
     Allocate (point(M,N), XX(M,N), Reg(Nlon,Nlat),Reg_pca(Nlon,Nlat))

     Do irow=1,N
        point(1:M,irow) = X(irow,1:M)
        XX(1:M,irow) = X(irow,1:M)
     Enddo
     Call Calc_all (Reg,N_Clus,"Clus_BL_"//trim(Dat_src)//"_",varn)
     Deallocate (point,XX)

    ! Calculate using pca of drought index:
    !Fname="/home/tanpv/Progs/Dr_SEA/EOF_New/pca_"//trim(Dat_src)//".nc"

    ! .True. if read in PCA, .False. if read in drought index
    !Call Read_data(Fname,"spe12",.true.)          ! Return: X(Nst,Ntime), Lon(Nst),Lat(Nst), Mask(Nlon,Nlat)
    !M=Ntime; N=Landpoint

     !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
     !Call Calc_all (Reg_pca,N_Clus,"Clus_pca_BL_"//trim(Dat_src)//"_",varn)

  !--------------------------------------------------------------------------
  !Do iclu=1,NClu_test        ! Test cases

  !Enddo

     
  CONTAINS
  !--------------------------------------------------------------------------
  Subroutine Calc_all (Reg,N_Clus,Fcase,varn)
     Character (Len=*) :: Fcase,varn
     Integer :: N_Clus, j,ist,ilon,ilat
     Real :: Reg(Nlon,Nlat)

     !N_Clus=Clus_test(iclu)  ! Num of clusters to be tested

     Allocate (Clus(Landpoint,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)
        Call test08 (point,M,N, cluster_num)
        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

     Call print_test
  
     ! For writing to file:
     Reg(:,:)=Rmiss   ! Regions
     ist=0
     Do ilat=1,Nlat
        Do ilon=1,Nlon
           if (Mask(ilon,ilat)/=Rmiss) then
              ist=ist+1
              Reg(ilon,ilat)=Clus(ist,N_Clus)
           endif
        Enddo
     Enddo
     !Call Write_Clus_to_NC (Reg,"Clus_BL_APH_",varn)
     Call Write_Clus_to_NC (Reg, trim(Fcase),varn)
     Deallocate (Clus,Clus_var,Clus_popu,Total_Ene,Total_var)
  !Enddo
  !Deallocate (X,Lon,Lat,Reg,xtime,point,XX)
  End  Subroutine Calc_all
  !stop 0
  !-----------------------------------------------------------------------------
  !CONTAINS
      !--------------------------------------------------------------------
      Subroutine Write_Clus_to_NC (Regs,Cases,varn)
        Character (Len=*) :: Cases,varn
        Character (Len=200) :: Fname
        Integer :: i,j
        Real :: tmp, tmp1, Percen_popu(N_Clus),Var(N_Clus)
        Real :: Regs(Nlon,Nlat)

        !write (Fname,"(A,I2.2,A)") "Clus_BL_APH_"//trim(varn)//"_",N_Clus,".nc"
        write (Fname,"(A,I2.2,A)") trim(Cases)//trim(varn)//"_",N_Clus,".nc"
        Call wr_ncdf_2d (Fname,xlon,xlat,nlon,nlat,Regs,"reg","Drought regions","-","lon","lat", Rmiss)
        tmp=Sum(Clus_popu(1:N_Clus,N_Clus))
        tmp1=Sum(Total_var(1:N_Clus))
        Do j=1,N_Clus
           Percen_popu(j)=Clus_popu(j,N_Clus)/tmp*100 
           Var(j)=Total_var(j)/tmp1*100  
        Enddo
        !write (Fname,"(A,I2.2,A)") "Clus_BL_APH_"//trim(varn)//"_",N_Clus,".txt"
        write (Fname,"(A,I2.2,A)") trim(Cases)//trim(varn)//"_",N_Clus,".txt"
        Open (3,File=Fname, Status="Unknown")
        write (3,*) N_Clus
        Do j=1,N_Clus
           write (3,"(2F6.2)") Percen_popu(j),Var(j)
        Enddo
        Close (3)
      End Subroutine Write_Clus_to_NC
      !--------------------------------------------------------------------
      Subroutine Read_data (Fname,varn,pca)
        Character (Len=*) :: Fname,varn
        Real, Allocatable :: PDSI(:,:,:), tmp(:), XX(:,:,:)
        Integer :: ist, ilon, ilat, ipoint, nt, it
        Logical :: pca

        !Fname="/home/tanpv/Progs/Dr_SEA/Results/BL_APH.nc"
        if (Allocated(xlon)) Deallocate (xlon)
        if (Allocated(xlat)) Deallocate (xlat)
        Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,nt,XX,trim(varn))

        Where (IsNaN(XX)); XX=Rmiss; Endwhere

        ! pca=true: For PCA of drought index:
        if (pca) then
           Ntime=Nt; it=1     ! Number of PCs
        else
           Ntime=Nt-Nmon; it=Nmon+1  ! Number of months
        endif
        Allocate (PDSI(Nlon,Nlat,Ntime))
        PDSI(:,:,1:Ntime) = XX(:,:,it:Nt)     ! Bo Nmon thang dau tien

        Nst=Nlon*Nlat
        if (.not.Allocated(Mask)) Allocate (Mask(Nlon,Nlat))
        Allocate (lon(Nst), lat(Nst), tmp(Ntime))
        ist=0; Landpoint=0; Mask(:,:)=Rmiss
        Do ilat=1, Nlat
           Do ilon=1, Nlon
              ist=ist+1;
              Lon(ist)=xlon(ilon); Lat(ist)=xlat(ilat)
              tmp(1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              if (.not.Ar1_Empty(tmp,Ntime)) then
                 Landpoint=Landpoint+1; Mask(ilon,ilat)=1
              endif
           Enddo
        Enddo
        if (.not.Allocated(X)) Allocate (X(Landpoint,Ntime))
        ipoint=0
        Do ilat=1, Nlat
           Do ilon=1, Nlon
              !tmp(1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              if (Mask(ilon,ilat)/=Rmiss) then
              !if (.not.Ar1_Empty(tmp,Ntime)) then
                 ipoint=ipoint+1
                 X(ipoint,1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              endif
           Enddo
        Enddo
        !print*,X; pause
        Deallocate (PDSI,XX,xtime,Lon,Lat,tmp)
      End Subroutine Read_data

      !--------------------------------------------------------------------
      Subroutine Read_data_xx (varn)
        Character (Len=200) :: Fname
        Real, Allocatable :: PDSI(:,:,:), tmp(:), XX(:,:,:)
        Integer :: ist, ilon, ilat, ipoint, nt, it
        Character (Len=*) :: varn

        !Fname="BL_VnGP.nc"   ! testing for Vietnam
        !Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,ntime,PDSI,"ori_pdsi")
        !Fname="pca_APH.nc"
        !Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,ntime,PDSI,"pdsi")

        Fname="/home/tanpv/Progs/Dr_SEA/Results/BL_APH.nc"
        !Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,ntime,PDSI,"ori_pdsi")
        Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,nt,XX,"spi_12")

        !Fname="/home/tanpv/Progs/Dr_SEA/EOF_New/pca_APH.nc"
        !Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,ntime,PDSI,"pdsi")

        !Fname="/home/tanpv/Progs/Dr_SEA/Results/BL_Ens.nc"
        !Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,ntime,PDSI,"ori_pdsi")

        ! For SPI, SPEI:
        if (trim(varn)=="spi_12".OR.trim(varn)=="spei_12") then
           Ntime=Nt-Nmon; it=Nmon+1
        else
           Ntime=Nt; it=1
        endif
        Allocate (PDSI(Nlon,Nlat,Ntime))
        PDSI(:,:,1:Ntime) = XX(:,:,it:Nt)     ! Bo Nmon thang dau tien

        Nst=Nlon*Nlat
        Allocate (Mask(Nlon,Nlat))
        Allocate (lon(Nst), lat(Nst), tmp(Ntime))
        ist=0; Landpoint=0; Mask(:,:)=Rmiss
        Do ilat=1, Nlat
           Do ilon=1, Nlon
              ist=ist+1;
              Lon(ist)=xlon(ilon); Lat(ist)=xlat(ilat)
              tmp(1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              if (.not.Ar1_Empty(tmp,Ntime)) then
                 Landpoint=Landpoint+1; Mask(ilon,ilat)=1
              endif
           Enddo
        Enddo
        Allocate (X(Landpoint,Ntime))
        ipoint=0
        Do ilat=1, Nlat
           Do ilon=1, Nlon
              !tmp(1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              if (Mask(ilon,ilat)/=Rmiss) then
              !if (.not.Ar1_Empty(tmp,Ntime)) then
                 ipoint=ipoint+1
                 X(ipoint,1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              endif
           Enddo
        Enddo
      End Subroutine Read_data_xx
   !------------------------------------------------------------------------------
   Logical Function Ar1_Empty(X,N)
     Integer :: N, i
     Real :: X(N)
     Logical :: Empty
     Empty=.false.
     Do i=1,N
        if (X(i)==Rmiss) then
           Empty=.true.; exit
        endif
     Enddo
     Ar1_Empty=Empty
   End Function Ar1_Empty
   !--------------------------------------------------------------------
  Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
  Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
  Include "kmean_lib.f90"
  !Include "Sub_kmean.f90"

!-----------------------------------------------------------------------------
 Subroutine print_test
 ! Clusters:
  !print*,"Clusters"
  !print"(3A7,20(3x,A1,I2.2))","STT","Lon","Lat",(("C",i),i=1,N_Clus)
  !Do i=1,Landpoint
  !   print"(I7,2F7.3,20I6)",i,Lon(i),Lat(i),(Clus(i,1:N_Clus))
  !Enddo

 !Populations:
  print*,"Populations"
  print"(A5,20(6x,A1,I2.2))","Gr",(("C",i),i=1,N_Clus)
  Do i=1,N_Clus
     print"(I5,20I9)", i, (Clus_popu(i,1:N_Clus))
  Enddo

 ! Variances:
  print*,"Variances"
  print"(A5,20(9x,A1,I2.2))","Gr",(("C",i),i=1,N_Clus)
  Do i=1,N_Clus
     print"(I5,20F12.1)", i, (Clus_var(i,1:N_Clus))
  Enddo

 ! Total variances:
  print*,"Mean variances"
  print"(A3,2A15)","Gr","Mean_Var" !,"Tot_Ene"
  Do i=1,N_Clus
     print"(I2,20F15.2)", i, Total_var(i) !, Total_Ene(i)
  Enddo
  End Subroutine print_test
!-----------------------------------------------------------------------------
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)

!*****************************************************************************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 )

  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)

!*****************************************************************************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

  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