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

#*****************************************************************************80
#
## python_mistake_test() displays some mistakes in python.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    24 October 2024
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'python_mistake_test():' )
  print ( '  python version: ' + platform.python_version ( ) )
  print ( '  numpy version:  ' + np.version.version )
  print ( '  Demonstrate some python mistakes.' )

  python_mistake01 ( )
  python_mistake02 ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'python_mistake_test():' )
  print ( '  Normal end of execution.' )

  return

def python_mistake01 ( ):

#*****************************************************************************80
#
## python_mistake01() tries to copy an array using the notation a = b.
#
#  Discussion:
#
#    Given arrays a and b, the statement "b = a" does not copy the entries
#    of a into b.  Instead, it allows b to be a second name for the information
#    in a.  A programmer who is unaware of this issue can set up disastrous
#    and mysterious errors.
#
#    By contrast, the correct procedure "b = a.copy()" does copy the entries
#    of a into b, creating two separate sets of memory.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    26 May 2019
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'python_mistake01():' )
  print ( '  Show that A=B is not the right way to copy arrays.' )

  print ( '' )
  print ( '  Use "=" to copy an array' )
  print ( '' )
  a = np.linspace ( 1, 10, 10 )
  print ( '  a =                        ', a )
  b = a
  print ( '  (Bad) b = a =              ', b )
  i1 = np.where ( b < 5 )
  b[i1] = -1
  i2 = np.where ( 5 <= b )
  b[i2] = +1
  print ( '  Do stuff to b, but do not touch a...')
  print ( '  a =                        ', a );
  print ( '  b =                        ', b );

  print ( '' )
  print ( '  Repeat, but use COPY rather than "="' )
  print ( '' )
  a = np.linspace ( 1, 10, 10 )
  print ( '  a =                        ', a )
  b = a.copy ( )
  print ( '  (Correct) b = a.copy() =   ', b )
  i1 = np.where ( b < 5 )
  b[i1] = -1
  i2 = np.where ( 5 <= b )
  b[i2] = +1
  print ( '  Do stuff to b, but do not touch a...')
  print ( '  a =                        ', a );
  print ( '  b =                        ', b );

  return

def python_mistake02 ( ):

#*****************************************************************************80
#
## python_mistake02(): Python assumes a data type for your array.
#
#  Discussion:
#
#    By specifying an array with integer values, Python may than assume
#    that integer arithmetic is to be applied.  This is disastrous when
#    division is involved.  In this example, the QR factors of an array
#    are to be computed, so that A = Q * R.  If A is specified using
#    integer values, a completely useless result is returned by the
#    function mgs_vector(), although the numpy function qr() is able to
#    process the data correctly.
#
#    The function mgs_vector() will work correctly if we simply specify
#    the array data using real number notation (decimal points!)
#
#  Modified:
#
#    24 October 2024
#
#  Author:
#
#    John Burkardt.
#
  import numpy as np

  print ( '' )
  print ( 'python_mistake02():' )
  print ( '  QR factorization on array stored as integer or real values.' )
  print ( '  User software fails for integer values.' )
#
#  Test 1:
#
  print ( '' )
  print ( '  A = np.array ( [ [1,2,3],[4,5,6],[7,8,10]]).' )
  print ( '  Algorithm     norm(Q^T*Q-I)    norm(Q*R-A)' )

  A = np.array ( [ [ 1, 2, 3 ], [ 4, 5, 6 ], [7,8,10]])
  m1, n1 = A.shape
  [ q3, r3 ] = mgs_vector ( A )
  [ q4, r4 ] = np.linalg.qr ( A )

  e1 = np.linalg.norm ( np.matmul ( q3.T, q3 ) - np.identity ( n1 ) )
  e2 = np.linalg.norm ( np.matmul ( q3, r3 ) - A, 'fro' )
  print ( '   mgs_vector    %e     %e' % ( e1, e2 ) )

  e1 = np.linalg.norm ( np.matmul ( q4.T, q4 ) - np.identity ( n1 ) )
  e2 = np.linalg.norm ( np.matmul ( q4, r4 ) - A, 'fro' )
  print ( '   np.linalg.qr  %e     %e' % ( e1, e2 ) )
#
#  Test 2:
#
  print ( '' )
  print ( '  A = np.array ( [ [1.0,2.0,3.0],[4.0,5.0,6.0],[7.0,8.0,10.0]]).' )
  print ( '  Algorithm     norm(Q^T*Q-I)    norm(Q*R-A)' )

  A = np.array ( [ [ 1.0, 2.0, 3.0 ], [ 4.0, 5.0, 6.0 ], [7.0,8.0,10.0]])
  m1, n1 = A.shape
  [ q3, r3 ] = mgs_vector ( A )
  [ q4, r4 ] = np.linalg.qr ( A )

  e1 = np.linalg.norm ( np.matmul ( q3.T, q3 ) - np.identity ( n1 ) )
  e2 = np.linalg.norm ( np.matmul ( q3, r3 ) - A, 'fro' )
  print ( '   mgs_vector    %e     %e' % ( e1, e2 ) )

  e1 = np.linalg.norm ( np.matmul ( q4.T, q4 ) - np.identity ( n1 ) )
  e2 = np.linalg.norm ( np.matmul ( q4, r4 ) - A, 'fro' )
  print ( '   np.linalg.qr  %e     %e' % ( e1, e2 ) )

  return

def mgs_vector ( A ):

#*****************************************************************************80
#
## mgs_vector() factors A=Q*R, using a vectorized algorithm.
#
#  Discussion:
#
#    Given m x n matrix A, the modified Gram-Schmidt algorithm computes
#    matrix factors Q and R such that:
#      A = Q*R, 
#    where Q is an m x min(m,n) orthogonal matrix
#    and R is an min(m,n) x n upper triangular matrix.
#
#  Author:
#
#    Dianne O'Leary, 09/2005
#
#  Reference: 
#
#    Golub and van Loan, Matrix Computations, Sec 5.2.8
#
  import numpy as np
#
#  Determine the matrix shape.
#
  m, n = A.shape
  Atype = A.dtype
#
#  Allocate storage.
#
  A2 = A.copy()
  A2 = A2.astype ( Atype )
  Q = np.zeros ( [ m, min ( m, n ) ], dtype = Atype )
  R = np.zeros ( [ min ( m, n ), n ], dtype = Atype )
#
#  Consider column k of A (k=1:min(m,n)),
#  Store it in column k of Q, after normalizing to length 1.
#  Orthogonalize each of the later columns of A against it,
#  saving the inner products in the k-th row of R.
#
  for k in range ( 0, min ( m, n ) ):
    R[k,k] = np.linalg.norm ( A2[:,k] )
    Q[:,k] = A2[:,k] / R[k,k]
    for j in range ( k + 1, n ):
      R[k,j] = np.dot ( np.conj ( Q[:,k] ), A2[:,j] )
      A2[:,j] = A2[:,j] - Q[:,k] * R[k,j]

  return Q, R

def timestamp ( ):

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

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

  return None

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