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

#*****************************************************************************80
#
## rref_cols_test() tests rref_cols().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    04 April 2025
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'rref_cols_test():' )
  print ( '  rref_cols() identifies the linearly independent columns' )
  print ( '  in a matrix, using the reduced row echelon form (RREF).' )

  A = np.array ( [ \
    [  1,  3,  0,  2,  6,  3,  1 ], \
    [ -2, -6,  0, -2, -8,  3,  1 ], \
    [  3,  9,  0,  0,  6,  6,  2 ], \
    [ -1, -3,  0,  1,  0,  9,  3 ] ] )
#
#  Use a wider output line so we can see the matrices properly.
#
  np.set_printoptions ( suppress = True, linewidth = 132 )

  print ( '' )
  print ( '  Matrix A:' )
  print ( A )

  A1, a_cols = rref_cols ( A )

  print ( '' )
  print ( '  rref_cols(A):' )
  print ( A1 )

  print ( '' )
  print ( '  Column indices:' )
  print ( a_cols )

  print ( '' )
  print ( '  Matrix columns:' )
  print ( A[:,a_cols] )

  return

def rref_cols ( A ):

#*****************************************************************************80
#
## rref_cols() identifies the independent columns in a matrix.
#
#  Discussion:
#
#    A rectangular matrix is in row reduced echelon form if:
#
#    * The leading nonzero entry in each row has the value 1.
#
#    * All entries are zero above and below the leading nonzero entry 
#      in each row.
#
#    * The leading nonzero in each row occurs in a column to
#      the right of the leading nonzero in the previous row.
#
#    * Rows which are entirely zero occur last.
#
#  Example:
#
#    M = 4, N = 7
#
#    Matrix A:
#
#     1    3    0    2    6    3    1
#    -2   -6    0   -2   -8    3    1
#     3    9    0    0    6    6    2
#    -1   -3    0    1    0    9    3
#
#    RREF(A):
#
#     1    3    0    0    2    0    0
#     0    0    0    1    2    0    0
#     0    0    0    0    0    1   1/3
#     0    0    0    0    0    0    0
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    27 December 2024
#
#  Author:
#
#    Original Python version by ChatGPT.
#    This version by John Burkardt.
#
#  Reference:
#
#    Charles Cullen,
#    An Introduction to Numerical Linear Algebra,
#    PWS Publishing Company, 1994,
#    ISBN: 978-0534936903,
#    LC: QA185.D37.C85.
#
#  Input:
#
#    real A(M,N), the matrix to be analyzed. 
#
#  Output:
#
#    real B(M,N), the reduced row echelon form of the matrix.
#
#    integer a_cols(*): the columns for which a pivot entry was found.
#
  import numpy as np
#
#  Initialize array of independent columns.
#
  a_cols = np.zeros ( 0, dtype = int )
#
#  Ensure the matrix elements are floating-point for division.
#
  A = A.astype ( float )
#
#  Get the array shape.
#
  rows, cols = A.shape
#
#  Set a tolerance for the size of a pivot value.
#
  tol = np.sqrt ( np.finfo(float).eps )
#
#  Start from the first row.
#
  row = 0
#
#  Seek a pivot for each column.
#
  for col in range ( cols ):
#
#  Exit if we have run out of rows to examine.
#
    if ( rows <= row ):
      break
#        
#  Find the pivot row.
#
    pivot_row = np.argmax ( np.abs ( A[row:rows, col] ) ) + row
#
#  Skip this column if all values are below the tolerance.
#
    if ( np.abs ( A[pivot_row, col] ) <= tol ):
      continue
#
#  A pivot was found for column col.
#
    a_cols = np.append ( a_cols, col )
#
#  Swap current row and pivot row.
#
    A[[row, pivot_row]] = A[[pivot_row, row]]
#  
#  Scale the pivot row to make the pivot element equal to 1.
#
    A[row] = A[row] / A[row, col]
#   
#  Eliminate the column below and above the pivot.
#
    for i in range ( rows ):
      if ( i != row ):
        A[i] = A[i] - A[i, col] * A[row]
#
#  Move to the next row.
#
    row = row + 1
    
  return A, a_cols

if ( __name__ == "__main__" ):
  rref_cols_test ( )