function svd_powers_test ( )

%*****************************************************************************80
%
%% svd_powers_test() tests svd_powers().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 December 2018
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../svd_powers' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'svd_powers_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test svd_powers().\n' );

  powers_data_test ( );
  svd_vectors_test ( );
  svd_project_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'svd_powers_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../svd_powers' )

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function powers_data_test ( )

%*****************************************************************************80
%
%% powers_data_test() tests powers_data().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 December 2018
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'powers_data_test():\n' );
  fprintf ( 1, '  powers_data() generates a dataset of powers.\n' );

  xmin = -1.0;
  xmax = +1.0;
  m = 11;
  n = 8;

  data = powers_data ( xmin, xmax, m, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  powers_data() generated this data:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : m
    for j = 1 : n
      fprintf ( 1, '  %10.4f', data(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end
function svd_project_test ( )

%*****************************************************************************80
%
%% svd_project_test() demonstrates projection using SVD vectors.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    06 December 2018
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'svd_project_test():\n' );
  fprintf ( 1, '  Compute singular vectors for x(i)^(j-1).\n' );
  fprintf ( 1, '  Then take an arbitrary vector.\n' );
  fprintf ( 1, '  Approximate it with 0, 1, 2, ... singular vectors.\n' );
%
%  Get the data array.
%
  xmin = -1.0;
  xmax = +1.0;
  m = 51;
  n = 8;

  x = linspace ( xmin, xmax, m );
  x = x';
  data = powers_data ( xmin, xmax, m, n );
%
%  Compute N SVD vectors.
%
  [ svd_vec, svd_val, avg ] = svd_vectors ( data, n );

  avg = avg / norm ( avg );
%
%  Get a data vector.
%
  se = sin ( x ) + 1 ./ ( x.^2 + 0.25 );
%
%  Break the data vector into:
%  SE0: constant + multiple of average
%  SE1: remainder.
%
  beta = ( se(m) - se(1) ) / ( avg(m) - avg(1) );
  alpha = se(1) - beta * avg(1);
  se0 = alpha + beta * avg;
  se1 = se - se0;
%
%  Construct approximations.
%
  fig_num = 0;

  for j = 0 : n

    fig_num = fig_num + 1;
    figure ( fig_num );
    clf ( );

    if ( j == 0 ) 
      sa = se0;
    else
      c = se1 * svd_vec(:,j)';
      sa = sa + c * svd_vec(:,j);
    end

    hold ( 'on' );
    plot ( x, se, 'b-', 'linewidth', 3 );
    plot ( x, sa, 'r-', 'linewidth', 3 );   
    grid ( 'on' );
    label = sprintf ( 'SVD projection using %d vectors\n', j );
    title ( label );
    hold ( 'off' );
    filename = sprintf ( 'svd_project_test_%d.png', j );
    print ( '-dpng', filename );
    fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  end

  return
end
function svd_vectors_test ( )

%*****************************************************************************80
%
%% svd_vectors_test() tests svd_vectors().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 December 2018
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'svd_vectors_test():\n' );
  fprintf ( 1, '  svd_vectors() computes SVD vectors.\n' );
  fprintf ( 1, '  Data is a set of power vectors.\n' );
  fprintf ( 1, '  Do SVD analysis.\n' );
  fprintf ( 1, '  Display the SVD vectors.\n' );
%
%  Get the data array.
%
  xmin = -1.0;
  xmax = +1.0;
  m = 51;
  n = 8;

  x = linspace ( xmin, xmax, m );
  data = powers_data ( xmin, xmax, m, n );
%
%  Compute N SVD vectors.
%
  [ svd_vec, svd_val, avg ] = svd_vectors ( data, n );
%
%  Display the average and the SVD vectors.
%
  fig_num = 0;

  for j = 0 : n
    fig_num = fig_num + 1;
    figure ( fig_num );
    clf
    if ( j == 0 )
      plot ( x, avg, 'Linewidth', 3 );
      label = 'SVD average';
    else
      plot ( x, svd_vec(:,j), 'Linewidth', 3 );
      label = sprintf ( 'SVD vector #%d  SVD value = %g', j, svd_val(j) );
    end
    title ( label, 'Fontsize', 16 );
    ylim ( [ -1.0, +1.0] );
    grid ( 'on' );
    filename = sprintf ( 'svd_vectors_%d.png', j );
    print ( '-dpng', filename );
    fprintf ( 1, '  Graphic saved as "%s"\n', filename );
  end
%
%  Display the singular values.
%
  fig_num = fig_num + 1;
  figure ( fig_num );
  clf ( );
  plot ( svd_val, 'ob-', 'Linewidth', 3 );
  title ( 'Singular values', 'Fontsize', 16 );
  grid ( 'on' );
  filename = 'svd_values.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphic saved as "%s"\n', filename );

  return
end

