subroutine cvtp_find_closest ( m, n, x, generator, width, modular, nearest )
!*****************************************************************************80
!
!! CVTP_FIND_CLOSEST finds the Voronoi cell generator closest to a point X.
!
! Discussion:
!
! This routine finds the closest Voronoi cell generator by checking every
! one. For problems with many cells, this process can take the bulk
! of the CPU time. Other approaches, which group the cell generators into
! bins, can run faster by a large factor.
!
! For this routine, if MODULAR is TRUE, then distance is done in a modular
! sense, as though the points were on a generalized torus. It's simple,
! really, we just need, in each coordinate, to consider
!
! X(I)-WIDTH(I), X(I), and X(I)+WIDTH(I).
!
! The bad part is, to keep our sanity, we want to replace X on output
! by the actual coordinates that got closest to some generator G,
! even though some of these coordinates may lie outside the unit
! hypercube. This is the right thing to do, so that the averaging
! process works correctly.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 25 July 2016
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer M, the spatial dimension.
!
! Input, integer N, the number of cell generators.
!
! Input, real ( kind = rk ) X(M), the point to be checked.
!
! Input, real ( kind = rk ) GENERATOR(M,N), the cell generators.
!
! Input, real ( kind = rk ) WIDTH(M), the width of the region in
! each dimension.
!
! Input, logical MODULAR, is TRUE if modular arithmetic is to be used.
!
! Output, integer NEAREST, the index of the nearest cell
! generators.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer m
integer n
real ( kind = rk ) generator(m,n)
real ( kind = rk ) dist_sq_min
real ( kind = rk ) dist_sq
integer i
integer j
logical modular
integer nearest
real ( kind = rk ) side
real ( kind = rk ) side1
real ( kind = rk ) side2
real ( kind = rk ) side3
real ( kind = rk ) width(m)
real ( kind = rk ) x(m)
real ( kind = rk ) y(m)
real ( kind = rk ) z(m)
nearest = 0
dist_sq_min = huge ( dist_sq_min )
do i = 1, n
dist_sq = 0.0D+00
do j = 1, m
if ( modular ) then
side1 = abs ( generator(j,i) - x(j) )
side2 = abs ( generator(j,i) + width(j) - x(j) )
side3 = abs ( generator(j,i) - width(j) - x(j) )
if ( side2 < side1 .and. side2 < side3 ) then
side = side2
y(j) = x(j) - width(j)
else if ( side3 < side1 .and. side3 < side2 ) then
side = side3
y(j) = x(j) + width(j)
else
side = side1
y(j) = x(j)
end if
else
side = abs ( generator(j,i) - x(j) )
y(j) = x(j)
end if
dist_sq = dist_sq + side ** 2
end do
if ( dist_sq < dist_sq_min ) then
dist_sq_min = dist_sq
nearest = i
z(1:m) = y(1:m)
end if
end do
!
! Overwrite X by Z, which is equal to X in modular arithmetic,
! but which is the closest to generator "NEAREST" (in non-modular
! arithmetic) of all the modularly equivalent copies of X.
!
x(1:m) = z(1:m)
return
end
subroutine cvtp_iteration ( m, n, generator, width, modular, &
sample_num_cvt, change_l2 )
!*****************************************************************************80
!
!! CVTP_ITERATION takes one step of the CVT iteration.
!
! Discussion:
!
! The routine is given a set of points, called "generators", which
! define a tessellation of the region into Voronoi cells. Each point
! defines a cell. Each cell, in turn, has a centroid, but it is
! unlikely that the centroid and the generator coincide.
!
! Each time this CVT iteration is carried out, an attempt is made
! to modify the generators in such a way that they are closer and
! closer to being the centroids of the Voronoi cells they generate.
!
! A large number of sample points are generated, and the nearest generator
! is determined. A count is kept of how many points were nearest to each
! generator. Once the sampling is completed, the location of all the
! generators is adjusted. This step should decrease the discrepancy
! between the generators and the centroids.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 05 Decemberc 2003
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer M, the spatial dimension.
!
! Input, integer N, the number of Voronoi cells.
!
! Input/output, real ( kind = rk ) GENERATOR(M,N), the Voronoi
! cell generators. On output, these have been modified
!
! Input, real ( kind = rk ) WIDTH(M), the width of the region in
! each direction.
!
! Input, logical MODULAR, is TRUE if modular arithmetic is to be used.
!
! Input, integer SAMPLE_NUM_CVT, the number of sample points.
!
! Output, real ( kind = rk ) CHANGE_L2, the sum of the L2 norms of the
! change in each generator's position.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer m
integer n
real ( kind = rk ) generator(m,n)
real ( kind = rk ) generator2(m,n)
real ( kind = rk ) change_gen
real ( kind = rk ) change_l2
integer count(n)
logical, parameter :: debug = .false.
integer i
integer j
logical modular
integer nearest
logical reset
integer sample_num_cvt
real ( kind = rk ) side
real ( kind = rk ) side1
real ( kind = rk ) side2
real ( kind = rk ) side3
real ( kind = rk ) width(m)
real ( kind = rk ) x(m)
generator2(1:m,1:n) = 0.0D+00
count(1:n) = 0
reset = .true.
do j = 1, sample_num_cvt
!
! Generate a sampling point X.
!
call cvtp_region_sampler ( m, 1, x, width )
reset = .false.
!
! Find the nearest cell generator G.
!
! Note that, to account for the modular arithemetic that is
! employed, the input value of X will be altered to
! the representative of X whose modular value is equal to X,
! but whose actual value is the closest to the generator G
! of all representatives of X.
!
! Otherwise, the averaging mechanism would be invalid!
!
call cvtp_find_closest ( m, n, x, generator, width, modular, nearest )
!
! Add X to the averaging data for GENERATOR(*,NEAREST).
!
generator2(1:m,nearest) = generator2(1:m,nearest) + x(1:m)
count(nearest) = count(nearest) + 1
end do
!
! Compute the new generators.
!
do j = 1, n
if ( count(j) /= 0 ) then
generator2(1:m,j) = generator2(1:m,j) / real ( count(j), kind = rk )
end if
end do
!
! It's possible that the generator would go outside the box.
! Use modular arithmetic to fix that.
!
if ( modular ) then
do j = 1, n
do i = 1, m
if ( generator2(i,j) < 0.0D+00 ) then
generator2(i,j) = generator2(i,j) + width(i)
else if ( width(i) < generator2(i,j) ) then
generator2(i,j) = generator2(i,j) - width(i)
end if
end do
end do
end if
!
! Determine the L2 norm of the change in the dataset.
!
! Because of our modular arithmetic, we need to do this carefully.
!
change_l2 = 0.0D+00
do j = 1, n
change_gen = 0.0D+00
if ( modular ) then
do i = 1, m
side1 = abs ( generator2(i,j) - generator(i,j) )
side2 = abs ( generator2(i,j) - width(i) - generator(i,j) )
side3 = abs ( generator2(i,j) + width(i) - generator(i,j) )
side = min ( side1, side2, side3 )
change_gen = change_gen + side ** 2
end do
else
do i = 1, m
side = abs ( generator2(i,j) - generator(i,j) )
change_gen = change_gen + side ** 2
end do
end if
change_l2 = change_l2 + sqrt ( change_gen )
end do
!
! Update.
!
generator(1:m,1:n) = generator2(1:m,1:n)
return
end
subroutine cvtp_region_sampler ( m, n, x, width )
!*****************************************************************************80
!
!! CVTP_REGION_SAMPLER returns a sample point in the physical region.
!
! Discussion:
!
! This code differs from the original CVT code only in that
! the WIDTH variable is available to specify the width of
! the box in each coordinate direction. Originally, the unit
! hypercube was used, and in fact, the current version of this
! program isn't really able to change the default widths from 1,
! but at least now, formally, the machinery is in place.
!
! This routine original interfaced with a lower routine called
! TEST_REGION, which tested whether the points generated in the
! bounding box were actually inside a possibly smaller physical
! region of interest. It's been a long time since that option
! was actually used, so it's been dropped.
!
! A point is chosen in the bounding box, either by a uniform random
! number generator, or from a vector Halton sequence.
!
! The original coding for this routine only supported a Halton
! sequence of dimension 3 or less. This restriction has been removed.
!
! Note that RESET was made an input-only quantity, in part to match
! the behavior of the routine in MATLAB, where it's cumbersome to
! support an input/output variable.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 05 December 2003
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer M, the spatial dimension.
!
! Input, integer N, the number of points to generate now.
!
! Output, real ( kind = rk ) X(M,N), the sample points.
!
! Input, real ( kind = rk ) WIDTH(M), the width of the region
! in each dimension.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer m
integer n
integer i
real ( kind = rk ) width(m)
real ( kind = rk ) x(m,n)
call random_number ( harvest = x(1:m,1:n) )
!
! Stretch the points to the given widths.
!
do i = 1, m
x(i,1:n) = width(i) * x(i,1:n)
end do
return
end
subroutine get_unit ( iunit )
!*****************************************************************************80
!
!! 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.
!
! If IUNIT = 0, then no free FORTRAN unit could be found, although
! all 99 units were checked (except for units 5, 6 and 9, which
! are commonly reserved for console I/O).
!
! 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.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 18 September 2005
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Output, integer IUNIT, the free unit number.
!
implicit none
integer i
integer ios
integer iunit
logical lopen
iunit = 0
do i = 1, 99
if ( i /= 5 .and. i /= 6 .and. i /= 9 ) 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 r8mat_write ( output_filename, m, n, table )
!*****************************************************************************80
!
!! R8MAT_WRITE writes an R8MAT file.
!
! Discussion:
!
! An R8MAT is an array of R8 values.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 31 May 2009
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) OUTPUT_FILENAME, the output file name.
!
! Input, integer M, the spatial dimension.
!
! Input, integer N, the number of points.
!
! Input, real ( kind = rk ) TABLE(M,N), the table data.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer m
integer n
integer j
character ( len = * ) output_filename
integer output_status
integer output_unit
character ( len = 30 ) string
real ( kind = rk ) table(m,n)
!
! Open the file.
!
call get_unit ( output_unit )
open ( unit = output_unit, file = output_filename, &
status = 'replace', iostat = output_status )
if ( output_status /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'R8MAT_WRITE - Fatal error!'
write ( *, '(a,i8)' ) ' Could not open the output file "' // &
trim ( output_filename ) // '" on unit ', output_unit
output_unit = -1
stop
end if
!
! Create a format string.
!
! For less precision in the output file, try:
!
! '(', m, 'g', 14, '.', 6, ')'
!
if ( 0 < m .and. 0 < n ) then
write ( string, '(a1,i8,a1,i8,a1,i8,a1)' ) '(', m, 'g', 24, '.', 16, ')'
!
! Write the data.
!
do j = 1, n
write ( output_unit, string ) table(1:m,j)
end do
end if
!
! Close the file.
!
close ( unit = output_unit )
return
end
subroutine timestamp ( )
!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
! Example:
!
! 31 May 2001 9:45:54.872 AM
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 18 May 2013
!
! Author:
!
! John Burkardt
!
implicit none
character ( len = 8 ) ampm
integer d
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
integer values(8)
integer y
call date_and_time ( values = 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 ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )
return
end