program main

!*****************************************************************************80
!
!! MAIN is the main program for CVT_MOD_DATASET.
!
!  Discussion:
!
!    CVT_MOD_DATASET computes a centroidal Voronoi Tessellation dataset, Mod 1.
!
!    This program is similar to the CVT_DATASET program, except that
!    the geometry of the region is not the unit square (or hypercube)
!    but rather the unit torus.  This is only in the sense of 
!    connectivity.  In other words, the top and bottom of the 
!    square behave as though they are "pasted together", as are 
!    the left and right sides.  In M dimensions, the effect is similar.
!    Really, all that is happening is that we are using modular arithmetic.
!
!    We hope, by this means, to be able to carry out CVT calculations
!    in which we do away with the boundary effects.  One sign of
!    boundary effects is the tendency of the CVT points to align
!    with the boundary, rather than to tend to the hexagonal arrangement
!    (in 2D at least).
!
!    It's nice to know we can do this, but we already know what the
!    points want to do in an infinite region in 2D, so is this just
!    a very expensive way to approximately compute points on a hexagonal
!    grid?  Well, yes, but it allows us to develop a tool that can
!    be used for other, related problems in which the solution is
!    NOT known.  In particular, for cases involving non-uniform densities
!    or various constraints, we might not know what the asymptotic
!    pattern of points is liable to be, and this code shows how to
!    design an algorithm that can avoid the effects of the boundary
!    lines, if they are, in fact, not an inherent part of the 
!    calculation.
!
!
!    The program is meant to be used interactively, or with
!    a simple script file.
!
!    Interaction is done through a series of keyword-driven commands.
!
!    The list of commands can be gotten by typing "HELP".
!
!    Briefly the user needs to specify the following:
!
!    * The spatial dimension M 
!      (this step may be omitted if the initial points are read from a file)
!
!      M = 2
!
!    * The number of points N
!      (this step may be omitted if the initial points are read from a file)
!
!      N = 100
!
!    * How the initial points are chosen:
!
!      INITIALIZE file_name
!        read the initial points from a file.
!        (this implicitly sets the spatial dimension M 
!        and the number of points N)
!      INITIALIZE RANDOM
!        initialize with N points in M dimensions, generated by RANDOM_NUMBER
!        (Fortran90 intrinsic);
!      INITIALIZE UNIFORM
!        initialize with N points in M dimensions, generated by UNIFORM;
!      INITIALIZE HALTON
!        initialize with N points in M dimensions, generated by HALTON;
!      INITIALIZE GRID
!        initialize with N points in M dimensions, generated by GRID;
!
!    * How the sampling is done:
!
!      SAMPLE_FUNCTION_CVT = RANDOM
!        generate sample points using RANDOM_NUMBER;
!      SAMPLE_FUNCTION_CVT = UNIFORM
!        generate sample points using UNIFORM;
!      SAMPLE_FUNCTION_CVT = HALTON
!        generate sample points using HALTON;
!      SAMPLE_FUNCTION_CVT = GRID
!        generate sample points using GRID;
!
!    * The number of sampling points to use:
!
!      SAMPLE_NUM_CVT = 200000
!
!    In the most common case, the user might have a set of data points
!    in a file, to be used as a starting point for the CVT iteration.
!
!    Currently, the program assumes that the region is contained in
!    the unit square.  
!
!  Example:
!
!    A reasonable set of input commands might be:
!
!      #  romero_01.inp
!      #
!      echo
!
!      m = 2
!      n = 100
!
!      initialize uniform
!      print
!
!      sample_function_cvt = uniform
!
!      sample_num_cvt = 5000
!      iterate 100
!      print
!
!      sample_num_cvt = 20000
!      iterate 5
!
!      sample_num_cvt = 80000
!      iterate 5
!
!      sample_num_cvt = 320000
!      iterate 5
!
!      write cvt.txt
!      quit
!
!    In the process shown here, the data file "data.txt" is read in.
!    The points in this file are used as starting points for the
!    cell generators.
!
!    The user then specifies various parameters, including that
!    sampling should be done with points generated by UNIFORM, that the current
!    generators should be printed out, that 100 iterations should
!    be made, with 5,000 sample points per generator.
!
!    After these 100 iterations, the sampling value is increased to
!    20,000, then 80,000, then 320,000 points, and each time,
!    5 more steps are taken.
!
!    Finally, the current generators are written to the file "cvt.txt".
!
!  Modified:
!
!    07 December 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    real ( kind = 8 ) CELL_GENERATOR(M,N), the Voronoi cell generators
!    of the Voronoi tessellation, as approximated by the CVT algorithm.  This
!    is the output quantity of most interest.
!
!    integer M, the spatial dimension.
!
!    integer MAXIT, the maximum number of correction iterations used in the
!    Voronoi calculation.  A typical value is 10.
!
!    integer N, the number of Voronoi cells to generate.
!
!    integer SAMPLE_FUNCTION_CVT, specifies how the region is sampled:
!    -1, the sampling function is RANDOM_NUMBER (Fortran90 intrinsic),
!    0, the sampling function is UNIFORM,
!    1, the sampling function is HALTON,
!    2, the sampling function is GRID.
!
!    integer SAMPLE_FUNCTION_INIT, specifies how the initial
!    generators are chosen:
!    -1, the initialization function is RANDOM_NUMBER (Fortran90 intrinsic),
!    0, the initialization function is UNIFORM,
!    1, the initialization function is HALTON,
!    2, the initialization function is GRID,
!    3, the initial values are read in from a file.
!
!    integer SAMPLE_NUM_CVT, the number of sampling points used on
!    each CVT iteration.  A typical value is 5000 * N.
!
!    integer SEED, determines how to initialize the random number routine.
!    If SEED is zero, then RANDOM_INITIALIZE will make up a seed
!    from the current real time clock reading.
!    If SEED is nonzero, then a reproducible sequence of random numbers
!    defined by SEED will be chosen.
!    By default, SEED initially has a value chosen by RANDOM_INITIALIZE,
!    but the user can reset SEED at any time.
!
  implicit none

  real ( kind = 8 ), allocatable, dimension ( :, : ) :: cell_generator
  real ( kind = 8 ) change_l2
  character ( len = 80 ) command
  character ( len = 80 ) command_cap
  integer m
  logical :: echo = .false.
  character ( len = 80 ) :: generator_file_in_name
  character ( len = 80 ) :: generator_file_out_name
  integer ierror
  integer it
  integer j
  integer last
  integer :: maxit = 10
  integer n
  logical reset
  integer :: sample_function_cvt = 0
  integer :: sample_function_init = 0
  integer sample_num_init
  integer :: sample_num_cvt = 10000
  logical s_eqi
  integer seed
  integer seed_init
  integer :: step = 0
  real ( kind = 8 ), allocatable, dimension ( : ) :: width
!
!  Print introduction and options.
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'CVT_MOD_DATASET'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Create Mod 1 CVT datasets.'
  write ( *, '(a)' ) '  Arithmetic: real (kind = 8 )'
  write ( *, '(a)' ) '  CVT = Centroidal Voronoi Tessellation.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The computational region is the unit hypercube.'
  write ( *, '(a)' ) '  Opposite boundaries are identified, yielding a'
  write ( *, '(a)' ) '  topology similar to a torus.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Given a region in ND, the problem is to determine'
  write ( *, '(a)' ) '  GENERATORS, a set of points which define a division'
  write ( *, '(a)' ) '  of the region into Voronoi cells, which are also'
  write ( *, '(a)' ) '  CENTROIDS of the Voronoi cells.'
  write ( *, '(a)' ) ' '
!
!  Get a default starting seed.
!
  call get_seed ( seed_init )

  seed = seed_init

  command = ' '

  do

    if ( command(1:1) == '#' ) then

    else
      write ( *, '(a,$)' ) ':>'
    end if

    read ( *, '(a)' ) command
    command_cap = command
    call s_cap ( command_cap )
    call s_blank_delete ( command_cap )

    if ( echo .or. command(1:1) == '#' ) then
      write ( *, '(a)' ) trim ( command )
    end if

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

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

    if ( command_cap(1:4) == 'ECHO' ) then

      echo = .true.
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'Command echo has been enabled.'

    else if ( command_cap(1:14) == 'INITIALIZEGRID' .or. &
            ( command_cap(1:15) == 'INITIALIZE=GRID' ) ) then

      sample_function_init = 2
      step = 0
      reset = .true.
      sample_num_init = n

      call region_sampler_mod ( m, n, sample_num_init, &
        cell_generator, sample_function_init, reset, seed, width )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The generators were initialized by GRID.'

    else if ( command_cap(1:16) == 'INITIALIZEHALTON' .or. &
              command_cap(1:17) == 'INITIALIZE=HALTON' ) then

      sample_function_init = 1
      step = 0
      reset = .true.
      sample_num_init = n
      
      call region_sampler_mod ( m, n, sample_num_init, &
        cell_generator, sample_function_init, reset, seed, width )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The generators were initialized by HALTON.'

    else if ( command_cap(1:16) == 'INITIALIZERANDOM' .or. &
              command_cap(1:17) == 'INITIALIZE=RANDOM' ) then

      sample_function_init = -1
      step = 0
      reset = .true.
      sample_num_init = n

      call region_sampler_mod ( m, n, sample_num_init, &
        cell_generator, sample_function_init, reset, seed, width )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The generators were initialized by RANDOM_NUMBER.'

    else if ( command_cap(1:17) == 'INITIALIZEUNIFORM' .or. &
              command_cap(1:18) == 'INITIALIZE=UNIFORM' ) then

      sample_function_init = 0
      step = 0
      reset = .true.
      sample_num_init = n

      call region_sampler_mod ( m, n, sample_num_init, &
        cell_generator, sample_function_init, reset, seed, width )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The generators were initialized by UNIFORM.'

    else if ( command_cap(1:1) == 'H' ) then

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'HELP!'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  M = value'
      write ( *, '(a)' ) '    Set the number of spatial dimensions.'
      write ( *, '(a)' ) '  N = value'
      write ( *, '(a)' ) '    Set the number of generators.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  ECHO'
      write ( *, '(a)' ) '    Turn command echo on.'
      write ( *, '(a)' ) '  HELP'
      write ( *, '(a)' ) '    Print this command list.'
      write ( *, '(a)' ) '  INITIALIZE GRID'
      write ( *, '(a)' ) '    Initializing function is GRID.'
      write ( *, '(a)' ) '  INITIALIZE HALTON'
      write ( *, '(a)' ) '    Initializing function is HALTON.'
      write ( *, '(a)' ) '  INITIALIZE RANDOM'
      write ( *, '(a)' ) '    Initializing function is RANDOM_NUMBER.'
      write ( *, '(a)' ) '  INITIALIZE UNIFORM'
      write ( *, '(a)' ) '    Initializing function is UNIFORM.'
      write ( *, '(a)' ) '  ITERATE n'
      write ( *, '(a)' ) '    Carry out N (more) iterations.'
      write ( *, '(a)' ) '  NOECHO'
      write ( *, '(a)' ) '    Cease command echo.'
      write ( *, '(a)' ) '  PRINT'
      write ( *, '(a)' ) '    Print problem parameters.'
      write ( *, '(a)' ) '  QUIT'
      write ( *, '(a)' ) '    Terminate execution.'
      write ( *, '(a)' ) '  READ = filename'
      write ( *, '(a)' ) '    Generators are initialized from a file.'
      write ( *, '(a)' ) '    (Automatically sets M and N.)'
      write ( *, '(a)' ) '  SAMPLE_FUNCTION_CVT = GRID'
      write ( *, '(a)' ) '    Sampling function is GRID.'
      write ( *, '(a)' ) '  SAMPLE_FUNCTION_CVT = HALTON'
      write ( *, '(a)' ) '    Sampling function is HALTON'
      write ( *, '(a)' ) '  SAMPLE_FUNCTION_CVT = RANDOM'
      write ( *, '(a)' ) '    Sampling function is RANDOM_NUMBER.'
      write ( *, '(a)' ) '  SAMPLE_FUNCTION_CVT = UNIFORM'
      write ( *, '(a)' ) '    Sampling function is UNIFORM.'
      write ( *, '(a)' ) '  SAMPLE_NUM_CVT = value'
      write ( *, '(a)' ) '    Set the number of CVT iteration sampling points.'
      write ( *, '(a)' ) '  SEED = value'
      write ( *, '(a)' ) '    Specify the random number seed.'
      write ( *, '(a)' ) '  WRITE filename'
      write ( *, '(a)' ) '    Generators are written from a file.'

    else if ( command_cap(1:7) == 'ITERATE' ) then

      if ( command_cap(8:8) == '=' ) then
        call s_to_i4 ( command_cap(9:), maxit, ierror, last )
      else if ( 0 < len_trim ( command_cap(8:) ) ) then
        call s_to_i4 ( command_cap(8:), maxit, ierror, last )
      else
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Enter MAXIT, the maximum number of'
        write ( *, '(a)' ) '  CVT iterations.'
        read ( *, * ) maxit
      end if

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  STEP  L2 Change'
      write ( *, '(a)' ) ' '

      do it = 1, maxit

        step = step + 1

        call cvt_iteration_mod ( m, n, cell_generator, width, sample_num_cvt, &
          sample_function_cvt, seed, change_l2 )

        write ( *, '(i4,f10.6)' ) step, change_l2

      end do

    else if ( command_cap(1:1) == 'M' ) then

      if ( command_cap(2:2) == '=' ) then
        call s_to_i4 ( command_cap(3:), m, ierror, last )
      else
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Enter M, the spatial dimension.'
        read ( *, * ) m
      end if

      if ( allocated ( cell_generator ) ) then
        deallocate ( cell_generator )
      end if

      if ( allocated ( width ) ) then
        deallocate ( width )
      end if

      allocate ( width(1:m) )
      width(1:m) = 1.0D+00

    else if ( command_cap(1:6) == 'NOECHO' ) then

      echo = .false.
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'Command echo has been disabled.'
!
!  Need to list "N" after "NOECHO" to avoid ambiguity.
!
    else if ( command_cap(1:1) == 'N' ) then

      if ( command_cap(2:2) == '=' ) then
        call s_to_i4 ( command_cap(3:), n, ierror, last )
      else
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Enter N, the number of generators.'
        read ( *, * ) n
      end if

      if ( allocated ( cell_generator ) ) then
        deallocate ( cell_generator )
      end if

      allocate ( cell_generator(1:m,1:n) )

    else if ( command_cap == 'PRINT' ) then

      call param_print ( m, n, maxit, seed, seed_init, &
        sample_function_cvt, sample_function_init, sample_num_cvt, &
        width )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'Current cell generators:'
      write ( *, '(a)' ) ' '
      do j = 1, n
        write ( *, * ) j, cell_generator(1:m,j)
      end do

    else if ( command_cap(1:1) == 'Q' ) then

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  QUIT requested.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CVT_MOD_DATASET'
      write ( *, '(a)' ) '  Normal end of execution.'
      exit

    else if ( s_eqi ( command(1:4), 'READ' ) ) then

      generator_file_in_name = command(5:)
      generator_file_in_name = adjustl ( generator_file_in_name )

      if ( generator_file_in_name(1:1) == '=' ) then
        generator_file_in_name(1:1) = ' '
        generator_file_in_name = adjustl ( generator_file_in_name )
      else if ( len_trim ( generator_file_in_name ) == 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Enter the generator file name.'
        read ( *, '(a)' ) generator_file_in_name
      end if

      sample_function_init = 3
      step = 0

      call file_column_count ( generator_file_in_name, m )

      call file_row_count ( generator_file_in_name, n )

      if ( allocated ( cell_generator ) ) then
        deallocate ( cell_generator )
      end if

      allocate ( cell_generator(1:m,1:n) )

      if ( allocated ( width ) ) then
        deallocate ( width )
      end if

      allocate ( width(1:m) )

      width(1:m) = 1.0D+00

      call data_read ( generator_file_in_name, m, n, cell_generator )

    else if ( command_cap(1:11) == 'SAMPLE_GRID' ) then

      sample_function_cvt = 2

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The region will be sampled using a fixed grid.'

    else if ( command_cap(1:19) == 'SAMPLE_FUNCTION_CVT' ) then

      if ( command_cap(20:20) /= '=' ) then

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Enter SAMPLE_FUNCTION_CVT:'
        write ( *, '(a)' ) '  -1: RANDOM_NUMBER (Fortran90 intrinsic);'
        write ( *, '(a)' ) '   0: UNIFORM;'
        write ( *, '(a)' ) '   1: HALTON;'
        write ( *, '(a)' ) '   2: GRID;'

        read ( *, * ) sample_function_cvt

      else

        if ( command_cap(21:26) == 'RANDOM' ) then
          sample_function_cvt = -1
        else if ( command_cap(21:27) == 'UNIFORM' ) then
          sample_function_cvt = 0
        else if ( command_cap(21:26) == 'HALTON' ) then
          sample_function_cvt = 1
        else if ( command_cap(21:24) == 'GRID' ) then
          sample_function_cvt = 2
        end if

      end if

      if ( sample_function_cvt == -1 ) then
        write ( *, '(a)' ) &
          '  Sampling function is RANDOM_NUMBER (Fortran90 intrinsic).'
      else if ( sample_function_cvt == 0 ) then
        write ( *, '(a)' ) '  Sampling function is UNIFORM.'
      else if ( sample_function_cvt == 1 ) then
        write ( *, '(a)' ) '  Sampling function is HALTON.'
      else if ( sample_function_cvt == 2 ) then
        write ( *, '(a)' ) '  Sampling function is GRID.'
      end if

    else if ( command_cap(1:14) == 'SAMPLE_NUM_CVT' ) then

      if ( command_cap(15:15) == '=' ) then
        call s_to_i4 ( command_cap(16:), sample_num_cvt, ierror, last )
      else
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) &
          '  Enter SAMPLE_NUM_CVT, the number of CVT sampling points.'
        read ( *, * ) sample_num_cvt
      end if

    else if ( command_cap(1:4) == 'SEED' ) then

      if ( command_cap(5:5) == '=' ) then
        call s_to_i4 ( command_cap(6:), seed, ierror, last )
      else
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Enter SEED, the initial random number seed.'
        read ( *, * ) seed
      end if

      write ( *, '(a,i12)' ) '  SEED = ', seed

      seed_init = seed

    else if ( command_cap(1:5) == 'WIDTH' ) then

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  I do not know how to read widths yet!'
      stop

    else if ( s_eqi ( command(1:5), 'WRITE' ) ) then

      generator_file_out_name = command(6:)
      generator_file_out_name = adjustl ( generator_file_out_name )

      if ( generator_file_out_name(1:1) == '=' ) then
        generator_file_out_name(1:1) = ' '
        generator_file_out_name = adjustl ( generator_file_out_name )
      else if ( len_trim ( generator_file_out_name ) == 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Enter the generator file name.'
        read ( *, '(a)' ) generator_file_out_name
      end if

      if ( len_trim ( generator_file_out_name ) == 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'Enter output file name:'
        read ( *, '(a)' ) generator_file_out_name
      end if

      call cvt_mod_write ( m, n, seed_init, sample_function_init, &
        generator_file_in_name, sample_function_cvt, sample_num_cvt, step, &
        change_l2, cell_generator, generator_file_out_name )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The generators were written to the file "' &
        // trim ( generator_file_out_name ) // '".'

    else

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Your input command was not recognized!'
      write ( *, '(a)' ) '  If ECHO is on, we will terminate.'

      if ( echo ) then
        exit
      end if

    end if

  end do

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end