program gene_cluster
!
!*******************************************************************************
!
!! GENE_CLUSTER groups gene expression data into clusters.
!
!
!  Discussion:
!
!    The data to be examined is assumed to be stored in a file.
!
!    The file is assumed to contain a number of records, with each
!    record stored on its own line.
!
!    Each record, in turn, contains a fixed number of data values
!    that describe a particular gene expression experiment.
!
!    Each record will be regarded as a point in N dimensional space.
!
!    The program will try to cluster the data, that is, to organize
!    the data by defining a number of cluster centers, which are
!    also points in N dimensional space, and assigning each record 
!    to the cluster associated with a particular center.
!
!    The method of assigning data aims to minimize the cluster energy,
!    which is taken to be the sum of the squares of the distances of
!    each data point from its cluster center.
!
!    In some contexts, it makes sense to use the usual Euclidean sort
!    of distance.  In others, it may make more sense to replace each
!    data record by a normalized version, and to assign distance
!    by computing angles between the unit vectors.
!
!  Modified:
!
!    12 September 2001
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer, parameter :: cluster_max = 10
  integer cluster_hi
  integer cluster_it_max
  integer cluster_lo
  integer cluster_num
  real, allocatable, dimension ( : ) :: cluster_size
  real, allocatable, dimension ( : ) :: cluster_size_inv
  character ( len = 100 ) data_file
  integer dim_num
  integer distance_type
  real, allocatable, dimension ( : ) :: energy
  integer energy_it_max
  integer i
  integer ierror
  integer inunit
  integer j
  integer normal
  real, allocatable, dimension ( :, :) :: point
  integer point_num
  integer seed
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'GENE_CLUSTER'
  write ( *, '(a)' ) '  Arrange gene expression data into clusters.'
  write ( *, '(a)' ) ' '

  seed = 1234567

  call random_initialize ( seed )
!
!  Get the name of the input data file.
!
  call s_input ( 'Enter the name of the data file', data_file, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GENE_CLUSTER'
    write ( *, '(a)' ) '  Error getting name of the data file.'
    write ( *, '(a)' ) '  Abnormal end of execution.'
    stop
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Data will be read from the file ' // trim ( data_file )
!
!  Determine the spatial dimension (columns) and number of points (rows).
!
  call file_column_count ( data_file, dim_num )
  call file_line_count ( data_file, point_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Spatial dimension is     ', dim_num
  write ( *, '(a,i6)' ) '  Number of data points is ', point_num
!
!  Allocate space for the POINT array.
!
  allocate ( point(1:dim_num,1:point_num) )
!
!  Read the data.
!
  open ( unit = inunit, file = data_file )
  do j = 1, point_num
    read ( inunit, * ) point(1:dim_num,j)
  end do
  close ( unit = inunit )
!
!  Get the range of cluster values to check.
!
  call i_range_input ( 'Enter lower and upper number of clusters', &
    cluster_lo, cluster_hi, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GENE_CLUSTER'
    write ( *, '(a)' ) '  Error getting cluster sizes.'
    write ( *, '(a)' ) '  Abnormal end of execution.'
    stop
  end if

  if ( cluster_hi > point_num ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  CLUSTER_HI exceeds number of points.'
    cluster_hi = point_num
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The program will try to determine the minimum energy'
  write ( *, '(a)' ) '  of a clustering, for cluster sizes in the range:'
  write ( *, '(2x,2i12)' ) cluster_lo, cluster_hi

  allocate ( cluster_size(1:cluster_hi) )
  allocate ( cluster_size_inv(1:cluster_hi) )
  allocate ( energy(1:cluster_hi) )

  do i = 1, cluster_hi
    cluster_size(i) = real ( i )
  end do
  cluster_size_inv(1:cluster_hi) = 1.0E+00 / cluster_size(1:cluster_hi)
!
!  Get the number of different random starts.
!
  call i_input ( 'Enter number of different random configurations to check', &
    cluster_it_max, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GENE_CLUSTER'
    write ( *, '(a)' ) '  Error getting number of random configurations.'
    write ( *, '(a)' ) '  Abnormal end of execution.'
    stop
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  For each number of clusters, the number of'
  write ( *, '(a)' ) '  distinct initial random configurations to be checked'
  write ( *, '(a,i6)' ) '  will be  ', cluster_it_max
!
!  Get the number of energy iterations for a particular random start:
!
  call i_input ( 'Enter the number of energy iterations', &
    energy_it_max, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'GENE_CLUSTER'
    write ( *, '(a)' ) '  Error getting number of energy iterations.'
    write ( *, '(a)' ) '  Abnormal end of execution.'
    stop
  end if
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  For each initial random configuration, the number of'
  write ( *, '(a)' ) '  times the program will recompute the cluster centers,'
  write ( *, '(a,i6)' ) '  cluster components, and energy is ', energy_it_max
!
!  Analyze the un-normalized data.
!
  normal = 0

  call analysis ( normal, dim_num, point_num, point, cluster_lo, &
    cluster_hi, cluster_it_max, energy_it_max, energy )

  cluster_num = cluster_hi + 1 - cluster_lo

  call data_to_gnuplot ( 'unnormal.txt', cluster_num, &
    cluster_size(cluster_lo:cluster_hi), energy(cluster_lo:cluster_hi) )

  call data_to_gnuplot ( 'unnormal2.txt', cluster_num, &
    cluster_size_inv(cluster_lo:cluster_hi), energy(cluster_lo:cluster_hi) )
!
!  Analyze the normalized data.
!
  normal = 1

  call analysis ( normal, dim_num, point_num, point, cluster_lo, & 
    cluster_hi, cluster_it_max, energy_it_max, energy )

  call data_to_gnuplot ( 'normal.txt', cluster_num, &
    cluster_size(cluster_lo:cluster_hi), energy(cluster_lo:cluster_hi) )

  call data_to_gnuplot ( 'normal2.txt', cluster_num, &
    cluster_size_inv(cluster_lo:cluster_hi), energy(cluster_lo:cluster_hi) )

  deallocate ( cluster_size )
  deallocate ( cluster_size_inv )
  deallocate ( energy )
  deallocate ( point )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'GENE_CLUSTER'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine analysis ( normal, dim_num, point_num, point, &
  cluster_lo, cluster_hi, cluster_it_max, energy_it_max, energy )
!
!*******************************************************************************
!
!! ANALYSIS computes the energy for a range of number of clusters.
!
!
!  Modified:
!
!    15 October 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NORMAL:
!    0, analyze the unnormalized data.
!    1, analyze the normalized data.
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of data points.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the data points.
!
!    Input, integer CLUSTER_LO, CLUSTER_HI, the low and high
!    cluster numbers that define a range of clusters to check.
!
!    Input, integer CLUSTER_IT_MAX, the number of different random
!    startup configurations to try.
!
!    Input, integer ENERGY_IT, the number of energy iterations.
!
  implicit none
!
  integer cluster_hi
  integer dim_num
  integer point_num
!
  real, dimension (dim_num,cluster_hi) :: cluster_center
  integer cluster_it
  integer cluster_it_max
  integer cluster_lo
  integer cluster_num
  real, dimension (cluster_hi) :: energy
  real energy_it
  real energy_max
  real energy_min
  integer energy_it_max
  integer i
  integer ierror
  integer j
  integer k
  real norm
  integer normal
  real, dimension (dim_num,point_num) :: point
  integer, dimension ( point_num ) :: point_cluster
  real, dimension (dim_num) :: r_min
  real, dimension (dim_num) :: r_max
  real t
!
!  If requested, normalize the data.
!
  if ( normal == 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Data is to be normalized.'
    write ( *, '(a)' ) ' '
    do j = 1, point_num
      norm = sqrt ( sum ( point(1:dim_num,j)**2 ) )
      if ( norm /= 0.0E+00 ) then
        point(1:dim_num,j) = point(1:dim_num,j) / norm
      end if
    end do
  end if

  if ( .false. ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  Point data:'
    write ( *, '(a)' ) ' '
    call point_print ( dim_num, point_num, point )
  end if
!
!  Compute the minimum and maximum component values.
!
  do i = 1, dim_num
    r_min(i) = minval ( point(i,1:point_num) )
    r_max(i) = maxval ( point(i,1:point_num) )
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Minimum and maximum data values:'
  write ( *, '(a)' ) ' '
  do i = 1, dim_num
    write ( *, '(i3,2f10.4)' ) i, r_min(i), r_max(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Cluster  Minimum    Maximum'
  write ( *, '(a)' ) '  Size     Energy     Energy'
  write ( *, '(a)' ) ' '
!
!  Consider a range of clusters.
!
  do cluster_num = cluster_lo, cluster_hi

    energy_min = huge ( energy_min )
    energy_max = 0.0E+00

    do cluster_it = 1, cluster_it_max

      call cluster_iteration ( dim_num, point_num, cluster_num, point, &
        r_min, r_max, cluster_center, energy_it, energy_it_max  )

      energy_min = min ( energy_min, energy_it )
      energy_max = max ( energy_max, energy_it )

    end do

    energy(cluster_num) = energy_min

    write ( *, '(i9,2f12.4)' ) cluster_num, energy_min, energy_max

  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Energy table:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Cluster              Energy'
  write ( *, '(a)' ) 'Size      Energy     /point    Sqrt(E/Pt)'
  write ( *, '(a)' ) ' '

  do cluster_num = cluster_lo, cluster_hi
    t = energy(cluster_num) / real ( point_num )
    write ( *, '(i9,3f12.4)' ) cluster_num, energy(cluster_num), t, sqrt ( t )
  end do

  return
end
subroutine cluster_iteration ( dim_num, point_num, cluster_num, point, &
  r_min, r_max, center, energy, energy_it_max )
!
!*******************************************************************************
!
!! CLUSTER_ITERATION seeks the minimal energy of a cluster of a given size.
!
!
!  Modified:
!
!    01 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of data points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the data points.
!
!    Input, real R_MIN(DIM_NUM), R_MAX(DIM_NUM), the coordinates of the
!    minimum and maximum corners of the region.
!
!    Output, real CENTER(DIM_NUM,CLUSTER_NUM), the centers associated with
!    the minimal energy clustering.
!
!    Output, real ENERGY, the total energy associated with the minimal 
!    energy clustering.
!
!    Input, integer ENERGY_IT, the number of energy iterations.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  real, dimension ( dim_num, cluster_num ) :: center
  real, dimension ( dim_num, cluster_num ) :: centroid
  integer, dimension ( point_num ) :: cluster
  real energy
  integer energy_it
  integer energy_it_max
  real energy_new
  integer i
  integer j
  integer missed
  real, dimension ( dim_num, point_num ) :: point
  real, dimension ( dim_num ) :: r
  real, dimension ( dim_num ) :: r_max
  real, dimension ( dim_num ) :: r_min
  integer, dimension ( cluster_num ) :: subs
!
!  Initialize the centers randomly.
!
  do j = 1, cluster_num
    call random_number ( harvest = r(1:dim_num) )
    center(1:dim_num,j) = ( 1.0E+00 - r(1:dim_num) ) * r_min(1:dim_num) &
                                    + r(1:dim_num)   * r_max(1:dim_num)
  end do
!
!  Initialize the clusters randomly.
!
  subs(1:cluster_num) = 0
  do i = 1, point_num
    call i_random ( 1, cluster_num, cluster(i) )
    j = cluster(i)
    subs(j) = subs(j) + 1
  end do

  call energy_computation ( dim_num, point_num, cluster_num, point, &
    center, cluster, energy )

  if ( .false. ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,g14.6)' ) &
      '  Initial total energy =                        ', energy
  end if

  if ( .false. ) then
    do j = 1, cluster_num
      write ( *, '(2i5,14f10.4)' ) j, subs(j), center(1:dim_num,j)
    end do
  end if

  do energy_it = 1, energy_it_max
!
!  #1:
!  Assign each point to the cluster of its nearest center.
!
    do i = 1, point_num
      call nearest_point ( dim_num, cluster_num, center, point(1,i), j )
      cluster(i) = j
    end do

    call energy_computation ( dim_num, point_num, cluster_num, point, &
      center, cluster, energy )

    if ( .false. ) then
      write ( *, '(a,g14.6)' ) &
        '  Total energy with old centers, new clusters = ', energy
    end if
!
!  #2:
!  Determine the centroids of the clusters.
!
    centroid(1:dim_num,1:cluster_num) = 0.0E+00
    subs(1:cluster_num) = 0

    do i = 1, point_num
      j = cluster(i)
      subs(j) = subs(j) + 1
      centroid(1:dim_num,j) = centroid(1:dim_num,j) + point(1:dim_num,i)
    end do

    missed = 0

    do j = 1, cluster_num

      if ( subs(j) /= 0 ) then
        centroid(1:dim_num,j) = centroid(1:dim_num,j) / real ( subs(j) )
      else
        missed = missed + 1

        centroid(1:dim_num,j) = point(1:dim_num,missed)

!       call random_number ( harvest = r(1:dim_num) )
!       centroid(1:dim_num,j) = ( 1.0E+00 - r(1:dim_num) ) * r_min(1:dim_num) &
!                                         + r(1:dim_num)   * r_max(1:dim_num)

      end if

    end do

    if ( missed /= 0 ) then
      if ( .false. ) then
        write ( *, '(a,i6)' ) '  Number of empty clusters was ', missed
      end if
    end if

    call energy_computation ( dim_num, point_num, cluster_num, point, &
      centroid, cluster, energy_new )

    if ( .false. ) then
      write ( *, '(a,g14.6)' ) &
        '  Total energy with new centers, new clusters = ', energy_new
    end if

    if ( .false. ) then
      do j = 1, cluster_num
        write ( *, '(2i6,14f10.4)' ) j, subs(j), center(1:dim_num,j)
        write ( *, '(12x,14f10.4)' )             centroid(1:dim_num,j)
      end do
    end if
!
!  Update the centers.
!
    center(1:dim_num,1:cluster_num) = centroid(1:dim_num,1:cluster_num)

  end do

  return
end
subroutine data_to_gnuplot ( file_name, n, x, y )
!
!*******************************************************************************
!
!! DATA_TO_GNUPLOT writes data to a file suitable for processing by GNUPLOT.
!
!
!  Discussion:
!
!    Once the data file is written, the GNUPLOT program can be used
!    to display the data, using commands like:
!
!      set term post default
!      set grid
!      set yrnage [0:*]
!      set title "Number of Clusters vs Total Energy, Normalized Data"
!      plot 'file_name'
!      quit
!
!  Modified:
!
!    12 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file
!    to which the data is written.
!
!    Input, integer N, the number of data values.
!
!    Input, real X(N), Y(N), the data.
!
  integer n
!
  character ( len = * ) file_name
  integer i
  integer iunit
  real x(n)
  real y(n)
!
  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'replace' )

  do i = 1, n
    write ( iunit, '(g16.8,2x,g16.8)' ) x(i), y(i)
  end do

  close ( unit = iunit )

  return
end
subroutine energy_computation ( dim_num, point_num, cluster_num, point, &
  center, cluster, energy )
!
!*******************************************************************************
!
!! ENERGY_COMPUTATION computes the total energy of a given clustering.
!
!
!  Discussion:
!
!    The total energy of a clustering is the sum of the energies of
!    each cluster.
!
!    The energy of a cluster is the sum of the squares of the distances of 
!    each point in the cluster to the center point of the cluster.
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of data points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the data points.
!
!    Input, real CENTER(DIM_NUM,CLUSTER_NUM), the center points.
!
!    Input, integer CLUSTER(POINT_NUM), indicates the cluster to which
!    each data point belongs.
!
!    Output, real ENERGY, the total energy of the clustering.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  real, dimension ( dim_num, cluster_num ) :: center
  integer, dimension ( point_num ) :: cluster
  real energy
  integer i
  integer j
  real, dimension ( dim_num, point_num ) :: point
!
  energy = 0.0E+00
  do i = 1, point_num
    j = cluster(i)
    energy = energy + sum ( ( center(1:dim_num,j) - point(1:dim_num,i) )**2 )
  end do

  return
end
subroutine file_column_count ( file_name, ncolumn )
!
!*******************************************************************************
!
!! FILE_COLUMN_COUNT counts the number of columns in the first line of a file.
!
!
!  Discussion:
!
!    The file is assumed to be a simple text file.
!
!    Most lines of the file is presumed to consist of NCOLUMN words, separated
!    by spaces.  There may also be some blank lines, and some comment lines,
!    which have a "#" in column 1.
!
!    The routine tries to find the first non-comment non-blank line and
!    counts the number of words in that line.
!
!    If all lines are blanks or comments, it goes back and tries to analyze
!    a comment line.
!
!  Modified:
!
!    21 June 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file.
!
!    Output, integer NCOLUMN, the number of columns assumed to be in the file.
!
  implicit none
!
  character ( len = * ) file_name
  logical got_one
  integer ios
  integer iunit
  character ( len = 256 ) line
  integer ncolumn
!
!  Open the file.
!
  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'old', form = 'formatted', &
    access = 'sequential', iostat = ios )

  if ( ios /= 0 ) then
    ncolumn = - 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Fatal error!'
    write ( *, '(a)' ) '  Could not open the file:'
    write ( *, '(a)' ) '    ' // trim ( file_name )
    return
  end if
!
!  Read one line, but skip blank lines and comment lines.
!
  got_one = .false.

  do

    read ( iunit, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if

    if ( len_trim ( line ) == 0 ) then
      cycle
    end if

    if ( line(1:1) == '#' ) then
      cycle
    end if

    got_one = .true.
    exit

  end do

  if ( .not. got_one ) then

    rewind ( iunit )

    do

      read ( iunit, '(a)', iostat = ios ) line

      if ( ios /= 0 ) then
        exit
      end if

      if ( len_trim ( line ) == 0 ) then
        cycle
      end if

      got_one = .true.
      exit

    end do

  end if

  close ( unit = iunit )

  if ( .not. got_one ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Warning!'
    write ( *, '(a)' ) '  The file does not seem to contain any data.'
    ncolumn = 0
    return
  end if

  call word_count ( line, ncolumn )

  return
end
subroutine file_line_count ( file_name, nline )
!
!*******************************************************************************
!
!! FILE_LINE_COUNT counts the number of lines in a file.
!
!
!  Discussion:
!
!    The file is assumed to be a simple text file.
!
!    Blank lines and comment lines, which begin with '#', are not counted.
!
!  Modified:
!
!    21 June 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file.
!
!    Output, integer NLINE, the number of lines found in the file.
!
  implicit none
!
  character ( len = * ) file_name
  integer ios
  integer iunit
  character ( len = 256 ) line
  integer nline
!
  nline = 0
!
!  Open the file.
!
  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'old', form = 'formatted', &
    access = 'sequential', iostat = ios )

  if ( ios /= 0 ) then
    nline = - 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILE_LINE_COUNT - Fatal error!'
    write ( *, '(a)' ) '  Could not open the file:' // trim ( file_name )
    return
  end if
!
!  Count the lines.
!
  do

    read ( iunit, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if

    if ( len_trim ( line ) == 0 ) then
      cycle
    end if

    if ( line(1:1) == '#' ) then
      cycle
    end if

    nline = nline + 1

  end do

  close ( unit = iunit )

  return
end
subroutine get_unit ( iunit )
!
!*******************************************************************************
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!
!  Discussion:
!
!    A "free" FORTRAN unit number is an integer between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!  Modified:
!
!    02 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5 and 6).
!
!    Otherwise, IUNIT is an integer between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
  implicit none
!
  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )

      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if

  end do

  return
end
subroutine i_input ( string, value, ierror )
!
!*******************************************************************************
!
!! I_INPUT prints a prompt string and reads an integer from the user.
!
!
!  Discussion:
!
!    If the input line starts with a comment character ('#') or is blank, 
!    the routine ignores that line, and tries to read the next one.
!
!  Modified:
!
!    27 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) STRING, the prompt string.
!
!    Output, integer VALUE, the value input by the user.
!
!    Output, integer IERROR, an error flag, which is zero if no error occurred.
!
  implicit none
!
  integer ierror
  integer last
  character ( len = 80 ) line
  character ( len = * ) string
  integer value
!
  ierror = 0
  value = huge ( value )
!
!  Write the prompt.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( string )

  do

    read ( *, '(a)', iostat = ierror ) line

    if ( ierror /= 0 ) then
      return
    end if
!
!  If the line begins with a comment character, go back and read the next line.
!
    if ( line(1:1) == '#' ) then
      cycle
    end if

    if ( len_trim ( line ) == 0 ) then
      cycle
    end if
!
!  Extract integer information from the string.
!
    call s_to_i ( line, value, ierror, last )

    if ( ierror /= 0 ) then
      value = huge ( value )
      return
    end if

    exit

  end do

  return
end
subroutine i_random ( ilo, ihi, i )
!
!*******************************************************************************
!
!! I_RANDOM returns a random integer in a given range.
!
!
!  Modified:
!
!    23 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ILO, IHI, the minimum and maximum acceptable values.
!
!    Output, integer I, the randomly chosen integer.
!
  implicit none
!
  integer i
  integer ihi
  integer ilo
  real r
  real, parameter :: rhi = 1.0E+00
  real, parameter :: rlo = 0.0E+00
!
  call r_random ( rlo, rhi, r )

  i = ilo + int ( r * real ( ihi + 1 - ilo ) )
!
!  In case of oddball events at the boundary, enforce the limits.
!
  i = max ( i, ilo )
  i = min ( i, ihi )

  return
end
subroutine i_range_input ( string, value1, value2, ierror )
!
!*******************************************************************************
!
!! I_RANGE_INPUT reads a pair of integers from the user, representing a range.
!
!
!  Discussion:
!
!    If the input line starts with a comment character ('#') or is blank,
!    the routine ignores that line, and tries to read the next one.
!
!  Modified:
!
!    27 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) STRING, the prompt string.
!
!    Output, integer VALUE1, VALUE2, the values entered by the user.
!
!    Output, integer IERROR, an error flag, which is zero if no error occurred.
!
  implicit none
!
  character, parameter :: comma = ','
  integer ierror
  integer last
  integer last2
  character ( len = 80 ) line
  character, parameter :: space = ' '
  character ( len = * ) string
  integer value1
  integer value2
!
  ierror = 0
  value1 = huge ( value1 )
  value2 = huge ( value2 )
!
!  Write the prompt.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( string )

  do

    read ( *, '(a)', iostat = ierror ) line

    if ( ierror /= 0 ) then
      return
    end if
!
!  If the line begins with a comment character, go back and read the next line.
!
    if ( line(1:1) == '#' ) then
      cycle
    end if

    if ( len_trim ( line ) == 0 ) then
      cycle
    end if
!
!  Remove commas.
!
    call s_rep_ch ( line, comma, space )
!
!  Extract integer information from the string.
!
    call s_to_i ( line, value1, ierror, last )

    if ( ierror /= 0 ) then
      value1 = huge ( value1 )
      return
    end if

    call s_to_i ( line(last+1:), value2, ierror, last2 )

    if ( ierror /= 0 ) then
      value2 = huge ( value2 )
      return
    end if

    exit

  end do

  return
end
subroutine nearest_point ( dim_num, cluster_num, center, point, nearest )
!
!*******************************************************************************
!
!! NEAREST_POINT finds the center point nearest a data point.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer CLUSTER_NUM, the number of center points.
!
!    Input, real CENTER(DIM_NUM,CLUSTER_NUM), the coordinates of the 
!    center points.
!
!    Input, real POINT(DIM_NUM), the data point to be checked.
!
!    Output, integer NEAREST, the index of the center point closest to
!    the data point.
!
  implicit none
!
  integer cluster_num
  integer dim_num
!
  real, dimension ( dim_num, cluster_num ) :: center
  real dist
  real dist_new
  integer j
  integer nearest
  real, dimension ( dim_num ) :: point
!
  dist = huge ( dist )
  nearest = 0

  do j = 1, cluster_num

    dist_new = sum ( ( point(1:dim_num) - center(1:dim_num,j) )**2 )

    if ( dist_new < dist ) then
      dist = dist_new
      nearest = j
    end if

  end do

  return
end
subroutine point_generate ( point_dist, dim_num, r_min, r_max, point_num, &
  point )
!
!*******************************************************************************
!
!! POINT_GENERATE generates data points for the problem.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer POINT_DIST, the point distribution to use.
!    1, use a uniform random distribution.
!    2, use a uniform grid of points.  (This hasn't been set up properly
!       except for 1D!).
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, real R_MIN(DIM_NUM), R_MAX(DIM_NUM), the coordinates of the
!    minimum and maximum corners of the region.
!
!    Input, integer POINT_NUM, the number of points to generate.
!
!    Output, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
  implicit none
!
  integer dim_num
  integer point_num
!
  integer i
  real point(dim_num,point_num)
  integer point_dist
  real r
  real r_max(dim_num)
  real r_min(dim_num)
!
  if ( point_dist == 1 ) then

    call random_number ( harvest = point(1:dim_num,1:point_num) )

    do i = 1, dim_num
      point(i,1:point_num) = r_min(i) + &
        point(i,1:point_num) * ( r_max(i) - r_min(i) )
    end do

  else if ( point_dist == 2 ) then

    do i = 1, point_num

      if ( point_num > 1 ) then
        r = real ( i - 1 ) / ( point_num - 1 )
      else
        r = 0.5E+00
      end if

      point(1:dim_num,i) = ( 1.0E+00 - r ) * r_min(1:dim_num) &
                                   + r   * r_max(1:dim_num)
    end do

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'POINT_GENERATE - Fatal error!'
    write ( *, '(a)' ) '  Meaningless input value of point distribution,'
    write ( *, '(a,i6)' ) '  POINT_DIST = ', point_dist
    stop

  end if

  return
end
subroutine point_print ( dim_num, point_num, point )
!
!*******************************************************************************
!
!! POINT_PRINT prints out the values of the data points.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
  implicit none
!
  integer dim_num
  integer point_num
!
  integer i
  real point(dim_num,point_num)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Data points:'
  write ( *, '(a)' ) ' '

  do i = 1, point_num
    write ( *, '(i6)' ) i
    write ( *, '(8f10.4)' ) point(1:dim_num,i)
  end do

  return
end
subroutine r_random ( rlo, rhi, r )
!
!*******************************************************************************
!
!! R_RANDOM returns a random real in a given range.
!
!
!  Modified:
!
!    06 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real RLO, RHI, the minimum and maximum values.
!
!    Output, real R, the randomly chosen value.
!
  implicit none
!
  real r
  real rhi
  real rlo
  real t
!
!  Pick T, a random number in (0,1).
!
  call random_number ( harvest = t )
!
!  Set R in ( RLO, RHI ).
!
  r = ( 1.0E+00 - t ) * rlo + t * rhi

  return
end
subroutine random_initialize ( seed )
!
!*******************************************************************************
!
!! RANDOM_INITIALIZE initializes the FORTRAN 90 random number seed.
!
!
!  Discussion:
!
!    If you don't initialize the random number generator, its behavior
!    is not specified.  If you initialize it simply by:
!
!      call random_seed
!
!    its behavior is not specified.  On the DEC ALPHA, if that's all you
!    do, the same random number sequence is returned.  In order to actually
!    try to scramble up the random number generator a bit, this routine
!    goes through the tedious process of getting the size of the random
!    number seed, making up values based on the current time, and setting
!    the random number seed.
!
!  Modified:
!
!    03 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer SEED.
!    If SEED is zero on input, then you're asking this routine to come up
!    with a seed value, which is returned as output.
!    If SEED is nonzero on input, then you're asking this routine to
!    use the input value of SEED to initialize the random number generator.
!
  implicit none
!
  integer count
  integer count_max
  integer count_rate
  integer i
  integer seed
  integer, allocatable :: seed_vector(:)
  integer seed_size
  real t
!
!  Initialize the random number seed.
!
  call random_seed
!
!  Determine the size of the random number seed.
!
  call random_seed ( size = seed_size )
!
!  Allocate a seed of the right size.
!
  allocate ( seed_vector(seed_size) )

  if ( seed /= 0 ) then

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'RANDOM_INITIALIZE'
    write ( *, '(a,i12)' ) '  Initialize RANDOM_NUMBER with user SEED = ', seed

  else

    call system_clock ( count, count_rate, count_max )

    seed = count

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'RANDOM_INITIALIZE'
    write ( *, '(a,i12)' ) &
      '  Initialize RANDOM_NUMBER with arbitrary SEED = ', seed

  end if
!
!  Now set the seed.
!
  seed_vector(1:seed_size) = seed

  call random_seed ( put = seed_vector(1:seed_size) )
!
!  Free up the seed space.
!
  deallocate ( seed_vector )
!
!  Call the random number routine a bunch of times.
!
  do i = 1, 100
    call random_number ( harvest = t )
  end do

  return
end
subroutine s_input ( string, value, ierror )
!
!*******************************************************************************
!
!! S_INPUT prints a prompt string and reads a string from the user.
!
!
!  Discussion:
!
!    If the input line starts with a comment character ('#') or is blank, 
!    the routine ignores that line, and tries to read the next one.
!
!  Modified:
!
!    27 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) STRING, the prompt string.
!
!    Output, character ( len = * ) VALUE, the value input by the user.
!
!    Output, integer IERROR, an error flag, which is zero if no error occurred.
!
  implicit none
!
  integer ierror
  character ( len = 80 ) line
  character ( len = * ) string
  character ( len = * ) value
!
  ierror = 0
!
!  Write the prompt.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( string )

  do

    read ( *, '(a)', iostat = ierror ) line

    if ( ierror /= 0 ) then
      value = 'S_INPUT: Input error!'
      return
    end if
!
!  If the line begins with a comment character, go back and read the next line.
!
    if ( line(1:1) == '#' ) then
      cycle
    end if

    if ( len_trim ( line ) == 0 ) then
      cycle
    end if

    value = line
    exit

  end do

  return
end
subroutine s_rep_ch ( s, c1, c2 )
!
!*******************************************************************************
!
!! S_REP_CH replaces all occurrences of one character by another.
!
!
!  Modified:
!
!    27 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character ( len = * ) S, the string.
!
!    Input, character C1, C2, the character to be replaced, and the
!    replacement character.
!
  implicit none
!
  character c1
  character c2
  integer i
  character ( len = * ) s
!
  do i = 1, len ( s )
    if ( s(i:i) == c1 ) then
      s(i:i) = c2
    end if
  end do

  return
end
subroutine s_to_i ( s, ival, ierror, last )
!
!*******************************************************************************
!
!! S_TO_I reads an integer value from a string.
!
!
!  Modified:
!
!    28 June 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, a string to be examined.
!
!    Output, integer IVAL, the integer value read from the string.
!    If the string is blank, then IVAL will be returned 0.
!
!    Output, integer IERROR, an error flag.
!    0, no error.
!    1, an error occurred.
!
!    Output, integer LAST, the last character of S used to make IVAL.
!
  implicit none
!
  character c
  integer i
  integer ierror
  integer isgn
  integer istate
  integer ival
  integer last
  character ( len = * ) s
!
  ierror = 0
  istate = 0
  isgn = 1
  ival = 0

  do i = 1, len_trim ( s )

    c = s(i:i)
!
!  Haven't read anything.
!
    if ( istate == 0 ) then

      if ( c == ' ' ) then

      else if ( c == '-' ) then
        istate = 1
        isgn = -1
      else if ( c == '+' ) then
        istate = 1
        isgn = + 1
      else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
        istate = 2
        ival = ichar ( c ) - ichar ( '0' )
      else
        ierror = 1
        return
      end if
!
!  Have read the sign, expecting digits.
!
    else if ( istate == 1 ) then

      if ( c == ' ' ) then

      else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
        istate = 2
        ival = ichar ( c ) - ichar ( '0' )
      else
        ierror = 1
        return
      end if
!
!  Have read at least one digit, expecting more.
!
    else if ( istate == 2 ) then

      if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
        ival = 10 * ival + ichar ( c ) - ichar ( '0' )
      else
        ival = isgn * ival
        last = i - 1
        return
      end if

    end if

  end do
!
!  If we read all the characters in the string, see if we're OK.
!
  if ( istate == 2 ) then
    ival = isgn * ival
    last = len_trim ( s )
  else
    ierror = 1
    last = 0
  end if

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end
subroutine word_count ( s, nword )
!
!*******************************************************************************
!
!! WORD_COUNT counts the number of "words" in a string.
!
!
!  Modified:
!
!    14 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, the string to be examined.
!
!    Output, integer NWORD, the number of "words" in the string.
!    Words are presumed to be separated by one or more blanks.
!
  implicit none
!
  logical blank
  integer i
  integer lens
  integer nword
  character ( len = * ) s
!
  nword = 0
  lens = len ( s )

  if ( lens <= 0 ) then
    return
  end if

  blank = .true.

  do i = 1, lens

    if ( s(i:i) == ' ' ) then
      blank = .true.
    else if ( blank ) then
      nword = nword + 1
      blank = .false.
    end if

  end do

  return
end