function p3_ex4 ( )

%*****************************************************************************80
%
%% p3_ex4, project 3, exercise 4
%
%  Discussion:
%
%    Clustering data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2019
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'p3_ex4:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Clustering of data.\n' );
%
%  Read the data.
%
  data = load ( 'faithful.txt' );
  x = data(:,1);
  y = data(:,2);
  n = length ( x );

  x_min = min ( x );
  x_max = max ( x );
  y_min = min ( y );
  y_max = max ( y );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of data values is %d\n', n );
  fprintf ( 1, '  %g <= x <= %g\n', x_min, x_max );
  fprintf ( 1, '  %g <= y <= %g\n', y_min, y_max );
%
%  Make a scatterplot of the (unnormalized) data.
%
  figure ( 1 )
  clf
  plot ( x, y, 'b.', 'Markersize', 15 );
  xlabel ( '<-- Duration (seconds) -->' )
  ylabel ( '<-- Wait (seconds) -->' );
  title ( 'Old Faithful eruption durations and waits', 'fontsize', 16 );
  grid ( 'on' );
  filename = 'p3_ex4_data.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Standardize the data.
%
  r = ( x - x_min ) / ( x_max - x_min );
  s = ( y - y_min ) / ( y_max - y_min );

  figure ( 2 )
  clf
  plot ( r, s, 'b.', 'Markersize', 15 );
  xlabel ( '<-- Duration (standardized) -->' )
  ylabel ( '<-- Wait (standarized) -->' );
  title ( 'Old Faithful eruption durations and waits', 'fontsize', 16 );
  grid ( 'on' );
  axis ( 'equal' );
  filename = 'p3_ex4_data_standardized.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Define two centers arbitrarily;
%  assign data to nearest.
%
  center = [ 0.40, 0.60; ...
             0.60, 0.30 ];
%
%  Determine which center is nearest.
%
  rs = [ r, s ];
  d1 = sum ( [ rs(:,1) - center(1,1), rs(:,2) - center(1,2) ].^2, 2 );

  rs1 = [ r - center(1,1), s - center(1,2) ];
  d1 = sum ( rs1.^2, 2 );
  rs2 = [ r - center(2,1), s - center(2,2) ];
  d2 = sum ( rs2.^2, 2 );
  c1 = find ( d1 < d2 );
  c2 = find ( d2 < d1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Cluster (%g,%g) has %d elements\n', center(1,1:2), length ( c1 ) );
  fprintf ( 1, '  Cluster (%g,%g) has %d elements\n', center(2,1:2), length ( c2 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Cluster energy = %g\n', sum ( d1(c1) ) );
  fprintf ( 1, '  Cluster energy = %g\n', sum ( d2(c2) ) );
  fprintf ( 1, '  Total energy   = %g\n',sum ( d1(c1) ) + sum ( d2(c2) ) );
%
%  Make a scatterplot of the clustered data.
%
  figure ( 3 )
  clf
  hold on
  plot ( r(c1), s(c1), 'b.', 'Markersize', 15 );
  plot ( center(1,1), center(1,2), 'c.', 'Markersize', 50 );
  plot ( r(c2), s(c2), 'r.', 'Markersize', 15 );
  plot ( center(2,1), center(2,2), 'm.', 'Markersize', 50 );
  xlabel ( '<-- Duration (standardized) -->' )
  ylabel ( '<-- Wait (standarized) -->' );
  title ( 'Old Faithful, initial clustering' );
  grid ( 'on' );
  axis ( 'equal' );
  hold off
  filename = 'p3_ex4_cluster_initial.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Update center by averaging.
%
  center(1,1) = mean ( r(c1) );
  center(1,2) = mean ( s(c1) );
  center(2,1) = mean ( r(c2) );
  center(2,2) = mean ( s(c2) );

  rs1 = [ r - center(1,1), s - center(1,2) ];
  d1 = sum ( rs1.^2, 2 );
  rs2 = [ r - center(2,1), s - center(2,2) ];
  d2 = sum ( rs2.^2, 2 );
  c1 = find ( d1 < d2 );
  c2 = find ( d2 < d1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Cluster (%g,%g) has %d elements\n', center(1,1:2), length ( c1 ) );
  fprintf ( 1, '  Cluster (%g,%g) has %d elements\n', center(2,1:2), length ( c2 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Cluster energy = %g\n', sum ( d1(c1) ) );
  fprintf ( 1, '  Cluster energy = %g\n', sum ( d2(c2) ) );
  fprintf ( 1, '  Total energy   = %g\n',sum ( d1(c1) ) + sum ( d2(c2) ) );
%
%  Make a scatterplot of the clustered data.
%
  figure ( 4 )
  clf
  hold on
  plot ( r(c1), s(c1), 'b.', 'Markersize', 15 );
  plot ( center(1,1), center(1,2), 'c.', 'Markersize', 50 );
  plot ( r(c2), s(c2), 'r.', 'Markersize', 15 );
  plot ( center(2,1), center(2,2), 'm.', 'Markersize', 50 );
  xlabel ( '<-- Duration (standardized) -->' )
  ylabel ( '<-- Wait (standarized) -->' );
  title ( 'Old Faithful, second clustering' );
  grid ( 'on' );
  axis ( 'equal' );
  hold off
  filename = 'p3_ex4_cluster_second.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Update center by kmeans.
%
  [ c, center ] = kmeans ( [ r, s ], 2 );
  c1 = c(1);
  c2 = c(2);

  rs1 = [ r - center(1,1), s - center(1,2) ];
  d1 = sum ( rs1.^2, 2 );
  rs2 = [ r - center(2,1), s - center(2,2) ];
  d2 = sum ( rs2.^2, 2 );
  c1 = find ( d1 < d2 );
  c2 = find ( d2 < d1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Cluster (%g,%g) has %d elements\n', center(1,1:2), length ( c1 ) );
  fprintf ( 1, '  Cluster (%g,%g) has %d elements\n', center(2,1:2), length ( c2 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Cluster energy = %g\n', sum ( d1(c1) ) );
  fprintf ( 1, '  Cluster energy = %g\n', sum ( d2(c2) ) );
  fprintf ( 1, '  Total energy   = %g\n',sum ( d1(c1) ) + sum ( d2(c2) ) );
%
%  Make a scatterplot of the clustered data.
%
  figure ( 5 )
  clf
  hold on
  plot ( r(c1), s(c1), 'b.', 'Markersize', 15 );
  plot ( center(1,1), center(1,2), 'c.', 'Markersize', 50 );
  plot ( r(c2), s(c2), 'r.', 'Markersize', 15 );
  plot ( center(2,1), center(2,2), 'm.', 'Markersize', 50 );
  xlabel ( '<-- Duration (standardized) -->' )
  ylabel ( '<-- Wait (standarized) -->' );
  title ( 'Old Faithful, k-means clustering' );
  grid ( 'on' );
  axis ( 'equal' );
  hold off
  filename = 'p3_ex4_cluster_kmeans.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'p3_ex4:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
