%
%  Warning: It is confusing, but the first column of X is actually the
%  2010 data, and the last column is 1890!
%
x = load ( 'snowfall.txt' );
%
%  Remove columns 1 and 10.
%
x = x(:,2:9);
%
%  Transpose.
%
x = x';
%
%  Compute SVD.
%
[u,s,v] = svd ( x );

if ( 0 )
  clf
  sdiag = diag ( s );
  plot ( 1:8, sdiag, 'b.-', 'Linewidth', 3, 'MarkerSize', 50 );
  title ( 'Singular values of SNOWFALL matrix' )
  grid on
end

if ( 0 )
  sdiag
  importance(1) = sdiag(1)^2 / sum ( sdiag.^2 )
  importance(2) = sdiag(2)^2 / sum ( sdiag.^2 )
  importance(3) = sdiag(3)^2 / sum ( sdiag.^2 )
end
%
%  Form the rank 1, 2, 3, 4, 5 approximants to A.
%
a1 = s(1,1) * u(:,1) * v(:,1)';
a2 = a1 + s(2,2) * u(:,2) * v(:,2)';
a3 = a2 + s(3,3) * u(:,3) * v(:,3)';
a4 = a3 + s(4,4) * u(:,4) * v(:,4)';
a5 = a4 + s(5,5) * u(:,5) * v(:,5)';

%
%  Column 1 of X is the 2010 snowfall.
%  Column 1 of A1 is the rank 1 approximant to 2010 snowfall.
%
if ( 0 )
clf
plot( t, x(:,1), 'r-', 'Linewidth', 3 )
grid on
title ( '2010 Snowfall' )
end

if ( 0 )
clf
hold on
plot( t, x(:,1), 'r-', 'Linewidth', 3 )
grid on
plot( t, a1(:,1), 'k-', 'Linewidth', 3 )
title ( 'Rank 1 Approx (black)' )
hold off

clf
hold on
plot( t, x(:,1), 'r-', 'Linewidth', 3 )
grid on
plot( t, a2(:,1), 'k-', 'Linewidth', 3 )
title ( 'Rank 2 Approx (black)' )
hold off

clf
hold on
plot( t, x(:,1), 'r-', 'Linewidth', 3 )
grid on
plot( t, a3(:,1), 'k-', 'Linewidth', 3 )
title ( 'Rank 3 Approx (black)' )
hold off

clf
hold on
plot( t, x(:,1), 'r-', 'Linewidth', 3 )
grid on
plot( t, a4(:,1), 'k-', 'Linewidth', 3 )
title ( 'Rank 4 Approx (black)' )
hold off

clf
hold on
plot( t, x(:,1), 'r-', 'Linewidth', 3 )
grid on
plot( t, a5(:,1), 'k-', 'Linewidth', 3 )
title ( 'Rank 5 Approx (black)' )
hold off

end
%
%  U contains the modes for variation over months in a fixed year.
%
if ( 0 )
us = r8col_normalize_li ( 8, 8, u );

clf
plot ( us(:,1), 'Linewidth', 3 )
grid on
title ( 'U Mode #1' )

plot ( us(:,2), 'Linewidth', 3 )
grid on
title ( 'U Mode #2' )

plot ( us(:,3), 'Linewidth', 3 )
grid on
title ( 'U Mode #3' )

plot ( us(:,4), 'Linewidth', 3)
grid on
title ( 'U Mode #4' )

end
%
%  V contains the modes for variation over years in a fixed month.
%
year = 1890:2010;
vs = r8col_normalize_li ( 121, 121, v );

%
%  Try to fix these plots by reversing the row ordering of v.
%
i = 121:-1:1;
j = 1:121;
vs(i,:) = vs(j,:);
if ( 1 )
  clf
  scatter ( year, vs(:,1), 'filled' )
  h = lsline
  set ( h, 'Linewidth', 3 )
  set ( h, 'Color', [0,0,0] )
  grid on
  xlabel ( '<--Year-->' )
  title ( 'V Mode #1' )
end

if ( 0 )
  clf
  scatter ( year, vs(:,2), 'filled' )
  h = lsline
  set ( h, 'Linewidth', 3 )
  set ( h, 'Color', [0,0,0] )
  grid on
  xlabel ( '<--Year-->' )
  title ( 'V Mode #2' )
end

if ( 0 )
  clf
  scatter ( year, vs(:,3), 'filled' )
  h = lsline
  set ( h, 'Linewidth', 3 )
  set ( h, 'Color', [0,0,0] )
  grid on
  xlabel ( '<--Year-->' )
  title ( 'V Mode #3' )
end

if ( 0 )
  clf
  scatter ( year, vs(:,4), 'filled' )
  h = lsline
  set ( h, 'Linewidth', 3 )
  set ( h, 'Color', [0,0,0] )
  grid on
  xlabel ( '<--Year-->' )
  title ( 'V Mode #4' )
end

