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

#*****************************************************************************80
#
## ex113, lab11 exercise 3.
#
#  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 ( 'ex113:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Read data about chemical properties of glass samples.' )
  print ( '  Copy data into matrix A.' )
  print ( '  Compute Amu = column averages.' )
  print ( '  Create centered matrix Ac = A - Amu.' )
  print ( '  Get SVD factorization of Ac.' )
  print ( '  Plot singular values S.' )
  print ( '  Plot cumulative sum of singular values.' )
  print ( '  Approximate A by .' )
#
#  Read data file.
#
  data = np.loadtxt ( 'glass_data.txt' )
#
#  Copy all but first and last columns of data into A.
#
  data_m, data_n = data.shape
  print ( '' )
  print ( '  data has ', data_m, ' rows and', data_n, 'columns.' )

  A = data[:,1:data_n-1]

# A = np.random.randn ( 6, 4 )

  m, n = A.shape
  print ( '' )
  print ( '  matrix A has ', m, ' rows and', n, 'columns.' )
#
#  Rewrite A = A0 + Ac
#
  mu = np.mean ( A, axis = 0 )
  A0 = np.zeros ( [ m, n ] )
  for j in range ( 0, n ):
    A0[0:m,j] = mu[j]

  Ac = np.zeros ( [ m, n ] )
  Ac = A - A0
#
#  Apply the SVD to Ac.
#
  U, svec, V = np.linalg.svd ( Ac, full_matrices = True )

  U_m, U_n = U.shape
  V_m, V_n = V.shape
#
#  Numpy made my life complicated by its peculiar SVD conventions!
#
  S_m = U_n
  S_n = V_m
  S = np.zeros ( [ S_n, S_m ] )
  S_num = min ( S_m, S_n )
  for i in range ( 0, S_num ):
    S[i,i] = svec[i]

  print ( '  matrix U has ', U_m, ' rows and', U_n, 'columns.' )
  print ( '  matrix S has ', S_m, ' rows and', S_n, 'columns.' )
  print ( '  matrix V has ', V_m, ' rows and', V_n, 'columns.' )
#
#  Plot the singular values.
#
  x = np.arange ( S_num )
  plt.plot ( x, svec, 'bo-', linewidth = 2 )
  plt.grid ( True )
  plt.xlabel ( '<-- I: index -->' )
  plt.ylabel ( '<-- Singular Values -->' )
  plt.title ( 'Singular values for glass data' )
  filename = 'ex113_1.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 ( S_num + 1 )
  y = [ 0.0 ]
  y = np.append ( y, np.cumsum ( svec**2 ) )
  y = y / np.sum ( svec**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 = 'ex113_2.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 = A0 + np.matmul ( U[:,0:1], np.matmul ( S[0:1,0:1], V[0:1,:] ) )
  e1 = np.linalg.norm ( A - A1 ) / A_norm
  print ( 'e1 = ', e1 )

  A2 = A0 + np.matmul ( U[:,0:2], np.matmul ( S[0:2,0:2], V[0:2,:] ) )
  e2 = np.linalg.norm ( A - A2 ) / A_norm
  print ( 'e2 = ', e2 )

  A3 = A0 + np.matmul ( U[:,0:3], np.matmul ( S[0:3,0:3], V[0:3,:] ) )
  e3 = np.linalg.norm ( A - A3 ) / A_norm
  print ( 'e3 = ', e3 )

  A4 = A0 + np.matmul ( U[:,0:4], np.matmul ( S[0:4,0:4], V[0:4,:] ) )
  e4 = np.linalg.norm ( A - A4 ) / A_norm
  print ( 'e4 = ', e4 )

  A5 = A0 + np.matmul ( U[:,0:5], np.matmul ( S[0:5,0:5], V[0:5,:] ) )
  e5 = np.linalg.norm ( A - A5 ) / A_norm
  print ( 'e5 = ', e5 )
#
#  Terminate.
#
  print ( '' )
  print ( 'ex113:' )
  print ( '  Normal end of execution.' )

  return

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