#! /usr/bin/env python3
#
def svd_circle ( A ):

#*****************************************************************************80
#
## svd_circle() plots the image of the unit circle under the 2x2 matrix A.
#
#  Discussion:
#
#    A typical 2x2 matrix will map the unit circle to some kind of tilted 
#    ellipse.  The aspect ratio of the ellipse (ratio of major to minor axes) 
#    is a measure of the conditioning of the matrix.
#
#    A singular matrix maps the unit circle to a line.
#
#    A diagonal matrix maps the unit circle to an ellipse with no tilting.
#    (There is no rotation.)
#
#    An orthogonal matrix maps the unit circle to the unit circle (but points
#    may have rotated.)
#
#    The identity matrix maps the unit circle to itself.
#
#    The singular value decomposition can predict the shape of the mapping:
#
#    * The vector V(1,1:2) maps to S(1)*U(1:2,1),
#    * The vector V(2,1:2) maps to S(2)*U(1:2,2).
#
#    The ellipse has an aspect ratio of S(1)/S(2).
#
#    The "tilt" or slope of the ellipse is U(2,1)/U(1,1).
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 March 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real A(2,2), the matrix whose mapping of the unit circle 
#    is to be studied.
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'svd_circle():' )
  print ( '  Given a matrix A,' )
  print ( '  plot points x on the unit circle, and' )
  print ( '  the images A*x of those points.' )
  print ( '  Show right singular vectors V' )
  print ( '  and their images, the scaled left singular vectors U.' )
  print ( '' )
  print ( '  Matrix A:' )
  print ( '' )
  print ( A )
#
#  This call produces matrices U, S and V.
#
  U, svec, V = np.linalg.svd ( A )

  print ( '  Singular values of A are', svec[0], 'and', svec[1] )
  print ( '  Right singular vectors are columns of V:' )
  V = V.T
  print ( V )
  print ( '  Left singular vectors are columns of U:' )
  print ( U )
  print ( '  Aspect ratio of ellipse will be', svec[0] / svec[1] )
  print ( '  Slope of ellipse will be', U[1,0] / U[0,0] )
#
#  Select N evenly spaced points X on the unit circle.
#
  n = 20

  i = np.linspace ( 0, n - 1, n, dtype = float )

  angle = 2.0 * np.pi * i / float ( n )

  r = 1.0

  x = np.array ( [ r * np.cos ( angle ), r * np.sin ( angle ) ] )
#
#  AX contains the image of each point under the mapping X -> A*X.
#
  ax = np.dot ( A, x )
#
#  Determine a common plot range.
#
  p_min = min ( np.min ( ax ), np.min ( x ) )
  p_max = max ( np.max ( ax ), np.max ( x ) )
#
#  Plot the corners of the region.  This is a trick to force MATLAB
#  to plot in a square with equal axes.  
#
  plt.clf ( )
  plt.scatter ( \
    [ p_min, p_max, p_max, p_min ], \
    [ p_min, p_min, p_max, p_max ] )
#
#  Label the points X on the unit circle.
#
  i = 97
  for j in range ( 0, n ):
    plt.text ( x[0,j], x[1,j], chr ( i ) )
    i = i + 1
    if ( 122 < i ):
      i = 97
#
#  Plot V1 and V2.
#
  plt.plot ( [ 0.0, V[0,0]], [ 0.0, V[1,0]], 'b', linewidth = 3 )
  plt.plot ( [ 0.0, V[0,1]], [ 0.0, V[1,1]], 'b', linewidth = 3 )

  plt.xlabel ( '--X axis--' )
  plt.ylabel ( '--Y axis--' )
  plt.title ( 'Points X on unit circle, and right singular vectors' )
  plt.axis ( np.array ( [ p_min, p_max, p_min, p_max ] ) )
  plt.axis ( 'equal' )
  plt.axis ( 'tight' )
  plt.axis ( 'square' )
  plt.grid ( True )
  filename = 'svd_circle_1.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )
#
#  Label the points AX on the image of the unit circle.
#
  plt.clf ( )
  plt.scatter ( np.array ( [ p_min, p_max, p_max, p_min ] ), \
                np.array ( [ p_min, p_min, p_max, p_max ] ) )
#
#  Plot the images AX of points X on the unit circle.
#
  i = 65
  for j in range ( 0, n ):
    plt.text ( ax[0,j], ax[1,j], chr ( i ) )
    i = i + 1
    if ( 90 < i ):
      i = 65
#
#  Plot the images of V1 and V2.
#
  plt.plot ( [ 0.0, svec[0]*U[0,0]], [ 0.0, svec[0]*U[1,0]], 'r', linewidth = 3 )
  plt.plot ( [ 0.0, svec[1]*U[0,1]], [ 0.0, svec[1]*U[1,1]], 'r', linewidth = 3 )

  plt.xlabel ( '--X axis--' )
  plt.ylabel ( '--Y axis--' )
  plt.title ( 'Images A*X of points on unit circle and scaled left singular vectors' )
  plt.axis ( np.array ( [ p_min, p_max, p_min, p_max ] ) )
  plt.axis ( 'equal' )
  plt.axis ( 'tight' )
  plt.axis ( 'square' )
  plt.grid ( True )
  filename = 'svd_circle_2.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'svd_circle():' )
  print ( '  Normal end of execution.' )

  return

def svd_circle_test ( ):

#*****************************************************************************80
#
## svd_circle_test() tests svd_circle().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 March 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'svd_circle_test():' )
  print ( '  Python version: ' + platform.python_version ( ) )
  print ( '  Test svd_circle().' )

  A = np.array ( [ \
    [ 3.0, 2.0 ],
    [ 1.0, 0.0 ] ] )

  svd_circle ( A )
#
#  Terminate.
#
  print ( '' )
  print ( 'svd_circle_test():' )
  print ( '  Normal end of execution.' )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  svd_circle_test ( )
  timestamp ( )