#! /usr/bin/env python3
#
def svd_glass ( ):

#*****************************************************************************80
#
## svd_glass
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    14 August 2019
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np
  import platform

  print ( '' )
  print ( 'svd_glass:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Read chemical properties of glass into matrix A.' )
  print ( '  Compute A0 = column averages.' )
  print ( '  Get SVD factorization of A - A0.' )
  print ( '  Plot singular values S.' )
  print ( '  Plot cumulative sum of singular values.' )
  print ( '  Approximate A by A0 + A1, A0 + A1 + A2, and so on' )
  print ( '  where A1 = S1 * U1 * V1, and so on.' )
#
#  Read the data file.
#
  A = np.loadtxt ( 'glass_data.txt' )
  m, n = A.shape
  print ( '' )
  print ( '  data has ', m, ' rows and', n, 'columns.' )
#
#  Compute A0, averages of each column of A, as an outer product.
#
  u0 = np.ones ( m )
  v0 = np.mean ( A, axis = 0 )
  A0 = np.outer ( u0, v0 )
#
#  Apply the SVD to A - A0.
#
  U, s, V = np.linalg.svd ( A - A0 )
#
#  Create a full matrix S from vector s.
#
  S = np.zeros ( A.shape )
  mn = min ( m, n )
  for i in range ( 0, mn ):
    S[i,i] = s[i]
#
#  Plot the singular values.
#
  x = np.arange ( mn )
  plt.plot ( x, s, 'bo-', linewidth = 2 )
  plt.grid ( True )
  plt.xlabel ( '<-- I: index -->' )
  plt.ylabel ( '<-- Singular Values -->' )
  plt.title ( 'Singular values for glass data' )
  filename = 'svd_glass_singular_values.png'
  plt.savefig ( filename )
  print ( '' )
  print ( '  Graphics saved in "%s"' % ( filename ) )
  plt.show ( )
  plt.clf ( )
#
#  Display the relative sum of the first K squared singular values.
#
  x = np.arange ( mn + 1 )
  y = [ 0.0 ]
  y = np.append ( y, np.cumsum ( s**2 ) )
  y = y / np.sum ( s**2 )
  plt.plot ( x, y, 'ro-', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<-- K: Number of components-->' )
  plt.ylabel ( '<-- Amount of variance explained-->' )
  plt.title ( 'Information proportion in first K components' )
  filename = 'svd_glass_percentage.png'
  plt.savefig ( filename )
  print ( '' )
  print ( '  Graphics saved in "%s"' % ( filename ) )
  plt.show ( )
  plt.clf ( )
#
#  A0 = mu
#
  A_norm = np.linalg.norm ( A )
  print ( '||A|| = ', A_norm )

  e0 = np.linalg.norm ( A - A0 ) / A_norm
  print ( 'e0 = ', e0 )

  A1 = s[0] * np.matmul ( U[:,0:1], V[0:1,:] )
  e1 = np.linalg.norm ( A - A0 - A1 ) / A_norm
  print ( 'e1 = ', e1 )

  A2 = s[1] * np.matmul ( U[:,1:2], V[1:2,:] )
  e2 = np.linalg.norm ( A - A0 - A1 - A2 ) / A_norm
  print ( 'e2 = ', e2 )

  A3 = s[2] * np.matmul ( U[:,2:3], V[2:3,:] )
  e3 = np.linalg.norm ( A - A0 - A1 - A2 - A3 ) / A_norm
  print ( 'e3 = ', e3 )

  A4 = s[3] * np.matmul ( U[:,3:4],  V[3:4,:] )
  e4 = np.linalg.norm ( A - A0 - A1 - A2 - A3 - A4 ) / A_norm
  print ( 'e4 = ', e4 )

  A5 = s[4] * np.matmul ( U[:,4:5], V[4:5,:] )
  e5 = np.linalg.norm ( A - A0 - A1 - A2 - A3 - A4 - A5 ) / A_norm
  print ( 'e5 = ', e5 )
#
#  Terminate.
#
  print ( '' )
  print ( 'svd_glass:' )
  print ( '  Normal end of execution.' )

  return

if ( __name__ == '__main__' ):
  svd_glass ( )