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

#*****************************************************************************80
#
## svd_6x4 uses SVD to carry out PCA on a random 6x4 matrix.
#
#  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_6x4:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Create a random 6x4 matrix.' )
  print ( '  Get the SVD factorization of A.' )

  m = 6
  n = 4
  A = np.random.randn ( 6, 4 )
  print ( '' )
  print ( '  Matrix A:' )
  print ( A )
#
#  Compute column averages:
#
  mu = np.mean ( A, axis = 0 )
  A0 = np.zeros ( A.shape )
  for j in range ( 0, n ):
    A0[0:m,j] = mu[j]
  Ac = A - A0
  print ( '' )
  print ( '  Norm of A =  ', np.linalg.norm ( A ) )
  print ( '  Norm of A0 = ', np.linalg.norm ( A0 ) )
  print ( '  Norm of Ac = ', np.linalg.norm ( Ac ) )
#
#  Get the singular value decomposition
#
  U, svec, V = np.linalg.svd ( Ac )
  S = np.zeros ( [ 6, 4] )
  for i in range ( 0, min ( 6, 4 ) ):
    S[i,i] = svec[i]
#
#  Verify the SVD factorization: A = A0 + U * S * V (no transpose!)
#
  USV = A0 + np.matmul ( U, np.matmul ( S, V ) )
  diff = np.linalg.norm ( A - USV )
  print ( '' )
  print ( '  Norm of A - U*S*V = %g' % ( diff ) )

  print ( '' )
#
#  Compute a 0 component approximation to A.
#
  diff = np.linalg.norm ( A - A0 )
  print ( '  Norm of A-A0 = %g' % ( diff ) )
#
#  Compute a 1 component approximation to A.
#
  A1 = A0 + np.matmul ( U[:,0:1], np.matmul ( S[0:1,0:1], V[0:1,:] ) )
  diff = np.linalg.norm ( A - A1 )
  print ( '  Norm of A-A1 = %g' % ( diff ) )
#
#  Compute a 2 component approximation to A.
#
  A2 = A0 + np.matmul ( U[:,0:2], np.matmul ( S[0:2,0:2], V[0:2,:] ) )
  diff = np.linalg.norm ( A - A2 )
  print ( '  Norm of A-A2 = %g' % ( diff ) )
#
#  Compute a 3 component approximation to A.
#
  A3 = A0 + np.matmul ( U[:,0:3], np.matmul ( S[0:3,0:3], V[0:3,:] ) )
  diff = np.linalg.norm ( A - A3 )
  print ( '  Norm of A-A3 = %g' % ( diff ) )
#
#  Compute a 4 component approximation to A.
#
  A4 = A0 + np.matmul ( U[:,0:4], np.matmul ( S[0:4,0:4], V[0:4,:] ) )
  diff = np.linalg.norm ( A - A4 )
  print ( '  Norm of A-A4 = %g' % ( diff ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'svd_6x4:' )
  print ( '  Normal end of execution.' )

  return

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