'''This module provides functions
re() and rre(), which compute the exact rational row echelon form and reduced row echelon form of an array of integers or rational numbers (fractions.Fraction),
printf() which prints an array of fractions (or integers or floats) nicely,
fractionized() which creates an array of rationals from an array of integers or integer-valued floats, and
decimalized() which creates an array of floats from an array of fractions or ints,
augmented() which is essentially hstack but with two arguments and the second may be a 1D array,
nullbasis() which provides an exact basis for the nullspace of an integer- or rational-valued matrix.
It also provides a class called "function" that supports arithmetic of numpy-friendly functions,
and "one" which is the constant function with value 1.
Dependencies: numpy, fractions, math (for gcd).
2/2/18 Added Matrix, M, etc to get rid of need for printf.
'''
'''
History ...
Formerly echelon.py.
Major reorganization 9/10/17. Now does all arithmetic with fractions.Fraction
and rejects any arrays with non-integer float elements.
I originally wrote this for datatype being *int* in particular,
then grafted on some code for handling Fraction.
'''
# 8/31/17
# Mistakes fixed. Some bullet-proofing added. I think it is correct.
# 8/30/17
# Add final step of converting element type to fractions.Fraction and make pivots 1
# JR 8/28/17
# generate exact rowechelon and redrowechelon for integer matrices
from numpy import *
from math import gcd
from fractions import Fraction
############### Functions for internal use not intended for use by end-users ###########################################################
def gcdl(l):
''' Not used by the latest re() and rre()
'''
g = l[0]
n = len(l)
i = 1
while i0 and a0.shape[1]>0:
a0type = type(a0[0,0])
if a0type is Fraction:
datatype = Fraction
elif a0type is float64:
datatype = float64
raise Exception('Row reduction of floats not reliable (zeros may not be exactly zero)')
m0,n0 = a1.shape
# The following is to create an array of python int rather than numpy.int64.
# The former permits arbitrarily large integers, while latter will wrap around, giving wrong answers.
A = array([[datatype(a1[i,j]) for j in range(n0)] for i in range(m0)],dtype=object) # an altered A will be returned by this function
a = A # a will be the current submatrix we're working on
while a.any(): # Note that any() on an array of *objects* gives weird (non-bool) answers.
# It seems to return the first non-zero value from the array or object zero.
# We are relying on any non-zero return value mapping to True, and 0 to False.
# Note that you can have an array with zero rows and/or zero columns. any() will return False for such arrays.
if a.shape[0] <= 1: break
j = argmax( array(a.any(axis=0),dtype=bool) > 0 ) # find the index of the leftmost non-zero column. There must be one.
a = a[:,j:] # shift right to omit all columns before j, which are already 0 and don't need to be touched. Column j becomes column 0.
if not ncols is None:
if a.shape[1] <= n0-ncols:
break
i = argmax( a[:,0] != 0 ) # find a pivot row. There must be one.
if i!=0: # otherwise no swap needed
a[[0,i]] = a[[i,0]] # swap the pivot row into row 0
m = a.shape[0] # number of rows in current submatrix
for i in range(1,m):
a[i] *= a[0,0] # a[0,0] is now necessarily non-zero
mu = a[i,0]//a[0,0]
a[i] -= mu*a[0]
# divide out gcd of each row
if a.shape[0]>0 and a.shape[1]>0 and type(a[0,0]) is int:
g = gcdl(a[i])
if g != 0:
a[i] //= gcdl(a[i])
if a[0,0] < 0: a[0] *= -1 # make the pivot element positive (This is a cleanup we can do within the integers.)
# we have now dealt with the top row, and have all zeros in first column below top row
a = a[1:,1:] # shift frame to deal with the next pivot
return A
def backred(a0,ncols=None):
'''This computes the reduced row echelon form of a 2D array of Fractions
that is already in integer row echelon form. It is used by rre() which is backred(re()).
Args: a0, a 2d numpy array of fractions.Fraction. Not altered by the function.
Returns: a new numpy array of fractions.Fraction.
'''
a = a0.copy()
m,n = a0.shape
for i in range(1,m)[::-1]:
#print('backred row',i)
if a[i].any():
j = argmax( a[i]!=0 ) # find pivot column. There must be one.
if not ncols is None and j >= ncols: continue
#print('\t pivot column',j, 'of',ncols,'considered')
for k in range(i):
#print('\t\tk =',k)
a[k] *= a[i,j]
mu = a[k,j]/a[i,j]
#print('zeroing in row',k)
a[k] -= mu*a[i]
# make all pivots 1
for i in range(m):
if a[i].any():
j = argmax( a[i]!=0 ) # find pivot column
if not ncols is None and j >= ncols: continue
a[i] /= a[i,j]
return a
############### The end-user functions ##############################################################################
def re(a0,ncols=None):
'''
This function reduces a 2D array of integers (any shape) to integer row echelon form.
Args:
a0 a 2D numpy array, either of numpy.int64 or int objects. Any shape is ok.
If a0 has any float elements, they will be rounded down to an integer. Perhaps an exception should be raised.
a0 is not altered by the function.
ncols, optional integer. If supplied, computation will cease once all pivots in the first ncols columns have been found.
Requested by Michael Cowen so that the when attempting to find inverse of a singular matrix by augmenting with the identity,
the answer agrees with that in the textbooks.
Returns:
a 2D numpy array of int objects, the integer row echelon form of a0. Not altered by the function.
'''
'''
This was derived from the original function now called rowechelonint().
'''
a1 = fractionized(a0)
m0,n0 = a1.shape
a = a1 # a will be the current submatrix we're working on, initially the whole thing
while a.any(): # Note that any() on an array of *objects* gives weird (non-bool) answers.
# It seems to return the first non-zero value from the array or object zero.
# We are relying on any non-zero return value mapping to True, and 0 to False.
# Note that you can have an array with zero rows and/or zero columns:
# any() will return False for such arrays.
if a.shape[0] == 0: break
j = argmax( array(a.any(axis=0),dtype=bool) > 0 ) # find the index of the leftmost non-zero column. There must be one.
a = a[:,j:] # shift right to omit all columns before j, which are already 0 and don't need to be touched. Column j becomes column 0.
if not ncols is None:
if a.shape[1] <= n0-ncols:
break
i = argmax( a[:,0] != 0 ) # find a pivot row. There must be one.
if i!=0: # otherwise no swap needed
a[[0,i]] = a[[i,0]] # swap the pivot row into row 0
m = a.shape[0] # number of rows in current submatrix
for i in range(1,m):
a[i] *= a[0,0] # a[0,0] is now necessarily non-zero
mu = a[i,0]/a[0,0]
a[i] -= mu*a[0]
a[0] /= a[0,0] # make the pivot 1
# we have now dealt with the top row, and have all zeros in first column below top row
a = a[1:,1:] # shift frame to deal with the next pivot
return Matrix(a1)
def det(a0,ncols=None):
# this is a modification of re() to keep track of row interchanges and row scalings
'''
This function reduces a 2D array of integers (any shape) to integer row echelon form.
Args:
a0 a 2D numpy array, either of numpy.int64 or int objects. Any shape is ok.
If a0 has any float elements, they will be rounded down to an integer. Perhaps an exception should be raised.
a0 is not altered by the function.
ncols, optional integer. If supplied, computation will cease once all pivots in the first ncols columns have been found.
Requested by Michael Cowen so that the when attempting to find inverse of a singular matrix by augmenting with the identity,
the answer agrees with that in the textbooks.
Returns:
a 2D numpy array of int objects, the integer row echelon form of a0. Not altered by the function.
'''
'''
This was derived from the original function now called rowechelonint().
'''
a1 = fractionized(a0)
m0,n0 = a1.shape
if m0!=n0:
print('Determinant not defined for a non-square matrix (',m0,n0,')')
return None
det = Fraction(1)
a = a1 # a will be the current submatrix we're working on, initially the whole thing
while a.any(): # Note that any() on an array of *objects* gives weird (non-bool) answers.
# It seems to return the first non-zero value from the array or object zero.
# We are relying on any non-zero return value mapping to True, and 0 to False.
# Note that you can have an array with zero rows and/or zero columns:
# any() will return False for such arrays.
if a.shape[0] == 0: break
j = argmax( array(a.any(axis=0),dtype=bool) > 0 ) # find the index of the leftmost non-zero column. There must be one.
if j>0: return Fraction(0)
a = a[:,j:] # shift right to omit all columns before j, which are already 0 and don't need to be touched. Column j becomes column 0.
if not ncols is None:
if a.shape[1] <= n0-ncols:
break
i = argmax( a[:,0] != 0 ) # find a pivot row. There must be one.
if i!=0: # otherwise no swap needed
a[[0,i]] = a[[i,0]] # swap the pivot row into row 0
det *= -1
#print('Swap')
m = a.shape[0] # number of rows in current submatrix
for i in range(1,m):
a[i] *= a[0,0] # a[0,0] is now necessarily non-zero
det /= a[0,0]
mu = a[i,0]/a[0,0]
a[i] -= mu*a[0]
det *= a[0,0]
a[0] /= a[0,0] # make the pivot 1
# we have now dealt with the top row, and have all zeros in first column below top row
a = a[1:,1:] # shift frame to deal with the next pivot
det *= a1.diagonal().prod()
if det.denominator==1:
return det.numerator
else:
return det
def rre(a,ncols=None):
'''Compute the reduced row echelon form of a 2D array of rational numbers.
Args: a, a 2d numpy array of fractions.Fraction, typically generated by re().
Not altered by the function.
Returns: a new numpy array of fractions.Fractions.
'''
return Matrix( backred(re(a,ncols),ncols) )
#from cfractions import complexFraction
#print(complexFraction(3,8))
def integerized(x):
'''
print(type(x))
print(x.__class__)
if isinstance(x,complexFraction): # this doesn't recognize object of type complexFraction
return x.integerized()
else:
print('hey')
return int(x)
'''
# Ugly temporary hack!
try:
return int(x)
except:
return x.integerized()
def fractionized(a): # might rename this "rationalized", or "fracked"?
'''This returns a new 2D numpy array of fraction.Fractions equivalent to the input array.
Args: a, a 2D numpy array of integers.
'''
a1 = array(a,dtype=object)
if len(a1.shape)==2:
# check that each element is exactly convertible to an int or is a Fraction
for i in range(a1.shape[0]):
for j in range(a1.shape[1]):
if not( type(a1[i,j]) is Fraction or integerized(a1[i,j])==a1[i,j] ):
raise Exception('fractionized() does not accept non-integer floating-point values such as '+str(a1[i,j]))
a1[i,j] = Fraction(a1[i,j])
else: # assuming 1D array
# check that each element is exactly convertible to an int or is a Fraction
for i in range(a1.shape[0]):
if not( type(a1[i]) is Fraction or integerized(a1[i])==a1[i] ):
raise Exception('fractionized() does not accept non-integer floating-point values such as '+str(a1[i]))
a1[i] = Fraction(a1[i])
return a1
def printf(a):
'''This pretty-prints a 2D numpy array of fraction.Fractions (or integers)
'''
if len(a.shape)==2:
stra = array( [[str(a[i,j]) for j in range(a.shape[1])] for i in range(a.shape[0])], dtype=object )
colw = array([max( [len(stra[i,j]) for i in range(a.shape[0])] ) for j in range(a.shape[1]) ])
colw += 1 # space
for i in range(a.shape[0]):
line = ''.join([ stra[i,j].rjust(colw[j]) for j in range(a.shape[1]) ])
print(line)
elif len(a.shape)==1:
print( ' '.join([ str(a[i]) for i in range(a.shape[0]) ] ) )
print()
def augmented(a,b):
'''
Return the augmented matrix a|b, where a is a 2D array and b is a 1D or 2D array.
The data type of the returned array will be inherited from a.
'''
#print(b.shape)
if len(b.shape)==1:
ab = empty( (a.shape[0],a.shape[1]+1), dtype=a.dtype )
ab[:,:-1] = a
ab[:,-1] = b
elif len(b.shape)==2:
assert( b.shape[0]==a.shape[0] )
ab = empty( (a.shape[0],a.shape[1]+b.shape[1]), dtype=a.dtype )
ab[:,:a.shape[1]] = a
ab[:,a.shape[1]:] = b
else:
raise Exception('augmented can\'t deal with this kind of b: '+str(b))
return Matrix(ab)
def decimalized(a):
return array( a, dtype=float )
from fractions import Fraction as fr
def nullbasis(a0): # Sep 14, 2017. This only works for a rational matrix
ra = rre( fractionized(a0) )
n = ra.shape[1]
pivotrows = where( ra.any(axis=1) )[0] # why does this return a tuple?
pivotcols = [argmax(ra[i]!=0) for i in pivotrows]
nonpivotcols = [j for j in range(n) if j not in pivotcols]
z = fractionized( zeros((n,n)) )
z[pivotcols] = -ra[pivotrows] # Note the negative sign!
fill_diagonal(z,fr(1))
b = z[:,nonpivotcols]
return Matrix(b)
class function:
'''A simple class for functional arithmetic.
'''
def __init__(self,f):
self.f = lambda x:f(x)
def __add__(self,other):
return function(lambda x: self.f(x) + other.f(x))
def __sub__(self,other):
return function(lambda x: self.f(x) - other.f(x))
def __mul__(self,other): # allows f*5
if isinstance(other,function): # Is this the only way to overload for multiple types??
return function(lambda x: self.f(x) * other.f(x))
else: # assume it's a scalar
return function(lambda x: other*self.f(x))
def __pow__(self,other): # allows f*5
if isinstance(other,function): # Is this the only way to overload for multiple types??
return function(lambda x: self.f(x) ** other.f(x))
else: # assume it's a scalar
return function(lambda x: self.f(x) ** other )
def __rmul__(self,other): # allows 5*f
if isinstance(other,function):
return function(lambda x: other.f(x) * self.f(x))
else: # assume it's a scalar
#print('scalar on left')
return function(lambda x: self.f(x)*other)
def __truediv__(self,other):
if isinstance(other,function):
return function(lambda x: self.f(x) / other.f(x))
else: # assume it's a scalar
return function(lambda x: self.f(x)/other )
def __call__(self,x):
return self.f(x)
def __abs__(self): # need for this not expected. quadrature requires it.
return function( lambda x: abs(self.f(x)) )
one = function(lambda x:1+0*x) # This might be nice to have.
#matrix = fractionized
#printm = printf
def sprintf(a):
'''This pretty-prints a 2D numpy array of fraction.Fractions (or integers)
'''
s = ''
if len(a.shape)==2:
stra = array( [[str(a[i,j]) for j in range(a.shape[1])] for i in range(a.shape[0])], dtype=object )
colw = array([max( [len(stra[i,j]) for i in range(a.shape[0])] ) for j in range(a.shape[1]) ])
colw += 1 # space
for i in range(a.shape[0]):
line = ''.join([ stra[i,j].rjust(colw[j]) for j in range(a.shape[1]) ])
s += line + '\n'
elif len(a.shape)==1:
s += ' '.join([ str(a[i]) for i in range(a.shape[0]) ] ) + '\n'
s += '\n'
return s
class M(ndarray):
def __repr__(self):
return sprintf(self)
def __str__(self):
return sprintf(self)
def Matrix(a):
return fractionized(a).view(M)
'''
if __name__ == '__main__':
a = random.randint(-3,9,(4,5))
print(a)
print()
printf(rre(a))
'''
'''
ringland@blue:~/public_html/309$ py echelon.py
[[6 1 2 0 6]
[2 5 1 4 2]
[6 4 6 3 3]
[5 7 2 1 4]]
1 0 0 0 589/453
0 1 0 0 -24/151
0 0 1 0 -124/151
0 0 0 1 115/453
ringland@blue:~/public_html/309$
'''