"""
PySCeS - Python Simulator for Cellular Systems (http://pysces.sourceforge.net)
Copyright (C) 2004-2009 B.G. Olivier, J.M. Rohwer, J.-H.S Hofmeyr all rights reserved,
Brett G. Olivier (bgoli@users.sourceforge.net)
Triple-J Group for Molecular Cell Physiology
Stellenbosch University, South Africa.
Permission to use, modify, and distribute this software is given under the
terms of the PySceS (BSD style) license. See LICENSE.txt that came with
this distribution for specifics.
NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK.
Brett G. Olivier
"""
from pysces.version import __version__
import time,os
import scipy
import scipy.linalg
import copy
__doc__ = """
PyscesStoich
------------
PySCeS stoichiometric analysis classes.
"""
## print 'Stoichiometry ver ' + __version__ + ' runtime: '+ time.strftime("%H:%M:%S")
__psyco_active__ = 0
## try:
## import psyco
## psyco.profile()
## __psyco_active__ = 1
## print 'PySCeS Stoichiometry Module is now PsycoActive!'
## except:
## __psyco_active__ = 0
##stoich_zero_valM = scipy.machar.MachAr().eps*2.0e4 # safer --> about 4e-11 (unofficially - a minpivot size)
##stoich_zero_valM = scipy.machar.MachAr().eps*2.0e4 # safe --> about 4e-12 (unofficially - a minpivot size)
##print '\tStoichiometric precision = ', stoich_zero_valM
class StructMatrix:
"""
This class is specifically designed to store structural matrix information
give it an array and row/col index permutations it can generate its own
row/col labels given the label src.
"""
array = None
ridx = None
cidx = None
row = None
col = None
shape = None
def __init__(self, array, ridx, cidx, row=None, col=None):
"""
Instantiate with array and matching row/col index arrays, optional label arrays
"""
self.array = array
self.ridx = ridx
self.cidx = cidx
self.row = row
self.col = col
self.shape = array.shape
def __call__(self):
return self.array
def getRowsByIdx(self, *args):
"""Return the rows referenced by index (1,3,5)"""
return self.array.take(args, axis=0)
def getColsByIdx(self, *args):
"""Return the columns referenced by index (1,3,5)"""
return self.array.take(args, axis=1)
def setRow(self, src):
"""
Assuming that the row index array is a permutation (full/subset)
of a source label array by supplying that source to setRow it
maps the row labels to ridx and creates self.row (row label list)
"""
self.row = [src[r] for r in self.ridx]
def setCol(self, src):
"""
Assuming that the col index array is a permutation (full/subset)
of a source label array by supplying that src to setCol
maps the row labels to cidx and creates self.col (col label list)
"""
self.col = [src[c] for c in self.cidx]
def getRowsByName(self, *args):
"""Return the rows referenced by label ('s','x','d')"""
assert self.row != None, "\nI need row labels"
try:
return self.array.take([self.row.index(l) for l in args], axis=0)
except Exception, ex:
print ex
print "\nValid row labels are: %s" % self.row
return None
def getColsByName(self, *args):
"""Return the columns referenced by label ('s','x','d')"""
assert self.col != None, "\nI need column labels"
try:
return self.array.take([self.col.index(l) for l in args], axis=1)
except Exception, ex:
print ex
print "Valid column labels are: %s" % self.col
return None
def getLabels(self, axis='all'):
"""Return the matrix labels ([rows],[cols]) where axis='row'/'col'/'all'"""
if axis == 'row': return self.row
elif axis == 'col': return self.col
else: return self.row, self.col
def getIndexes(self, axis='all'):
"""Return the matrix indexes ([rows],[cols]) where axis='row'/'col'/'all'"""
if axis == 'row': return self.ridx
elif axis == 'col': return self.cidx
else: return self.ridx, self.cidx
def getByIdx(self, row, col):
assert row in self.ridx, '\n%s is an invalid index' % row
assert col in self.cidx, '\n%s is an invalid index' % col
return self.array[row, col]
def getByName(self, row, col):
assert row in self.row, '\n%s is an invalid name' % row
assert col in self.col, '\n%s is an invalid name' % col
return self.array[self.row.index(row), self.col.index(col)]
def setByIdx(self, row, col, val):
assert row in self.ridx, '\n%s is an invalid index' % row
assert col in self.cidx, '\n%s is an invalid index' % col
self.array[row, col] = val
def setByName(self, row, col, val):
assert row in self.row, '\n%s is an invalid name' % row
assert col in self.col, '\n%s is an invalid name' % col
self.array[self.row.index(row), self.col.index(col)] = val
class MathArrayFunc(object):
"""A class of basic array functions some LAPACK based"""
__doc__ = '''PySCeS array functions - used by Stoich'''
array_kind = {'i':0, 'l': 0, 'f': 0, 'd': 0, 'F': 1, 'D': 1}
array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
array_type = [['f', 'd'], ['F', 'D']]
LinAlgError = 'LinearAlgebraError'
def commonType(self,*arrays):
"""
commonType(\*arrays)
Numeric detect and set array precision (will be replaced with new scipy.core compatible code when ready)
Arguments:
\*arrays: input arrays
"""
kind = 0
# precision = 0
# force higher precision in lite version
precision = 1
for a in arrays:
t = a.dtype.char
kind = max(kind, self.array_kind[t])
precision = max(precision, self.array_precision[t])
return self.array_type[kind][precision]
def castCopyAndTranspose(self,type, *arrays):
"""
castCopyAndTranspose(type, \*arrays)
Cast numeric arrays to required type and transpose
Arguments:
type: the required type to cast to
\*arrays: the arrays to be processed
"""
cast_arrays = ()
for a in arrays:
if a.typecode() == type:
cast_arrays = cast_arrays + (copy.copy(scipy.transpose(a)),)
else:
cast_arrays = cast_arrays + (copy.copy(scipy.transpose(a).astype(type)),)
if len(cast_arrays) == 1:
return cast_arrays[0]
else:
return cast_arrays
def assertRank2(self,*arrays):
"""
assertRank2(\*arrays)
Check that we are using a 2D array
Arguments:
\*arrays: input array(s)
"""
for a in arrays:
if len(a.shape) != 2:
raise LinAlgError, 'Array must be two-dimensional'
def SwapCol(self,res_a,r1,r2):
"""
SwapCol(res_a,r1,r2)
Swap two columns using BLAS swap, arrays can be (or are upcast to) type double (d) or double complex (D).
Returns the colswapped array
Arguments:
res_a: the input array
r1: the first column to be swapped
r2: the second column to be swapped
"""
if self.array_kind[self.commonType(res_a)] == 1: # brett 20041226 added complex support to swap
res_a = res_a.astype('D')
return self.SwapColz(res_a,r1,r2)
else:
res_a = res_a.astype('d')
return self.SwapCold(res_a,r1,r2)
def SwapRow(self,res_a,r1,r2):
"""
SwapRow(res_a,r1,r2)
Swaps two rows using BLAS swap, arrays can be (or are upcast to) type double (d) or double complex (D).
Returns the rowswapped array.
Arguments:
res_a: the input array
r1: the first row index to be swapped
r2: the second row index to be swapped
"""
if self.array_kind[self.commonType(res_a)] == 1: # brett 20041226 added complex support to swap
res_a = res_a.astype('D')
return self.SwapRowz(res_a,r1,r2)
else:
res_a = res_a.astype('d')
return self.SwapRowd(res_a,r1,r2)
def SwapElem(self,res_a,r1,r2):
"""
SwapElem(res_a,r1,r2)
Swaps two elements in a 1D vector
Arguments:
res_a: the input vector
r1: index 1
r2: index 2
"""
res_a[r1],res_a[r2] = res_a[r2],res_a[r1]
return (res_a)
def SwapCold(self,res_a,c1,c2):
"""
SwapCold(res_a,c1,c2)
Swaps two double (d) columns in an array using BLAS DSWAP. Returns the colswapped array.
Arguments:
res_a: input array
c1: column index 1
c2: column index 2
"""
res_a = res_a.astype('d')
res_a[:,c1],res_a[:,c2] = scipy.linalg.fblas.dswap(res_a[:,c1],res_a[:,c2])
return (res_a)
def SwapRowd(self,res_a,r1,r2):
"""
SwapRowd(res_a,c1,c2)
Swaps two double (d) rows in an array using BLAS DSWAP. Returns the rowswapped array.
Arguments:
res_a: input array
c1: row index 1
c2: row index 2
"""
res_a = res_a.astype('d')
res_a[r1,:],res_a[r2,:] = scipy.linalg.fblas.dswap(res_a[r1,:],res_a[r2,:])
return (res_a)
def SwapColz(self,res_a,c1,c2):
"""
SwapColz(res_a,c1,c2)
Swaps two double complex (D) columns in an array using BLAS ZSWAP. Returns the colswapped array.
Arguments:
res_a: input array
c1: column index 1
c2: column index 2
"""
res_a = res_a.astype('D')
res_a[:,c1],res_a[:,c2] = scipy.linalg.fblas.zswap(res_a[:,c1],res_a[:,c2])
return (res_a)
def SwapRowz(self,res_a,r1,r2):
"""
SwapRowz(res_a,c1,c2)
Swaps two double complex (D) rows in an array using BLAS ZSWAP. Returns the rowswapped array.
Arguments:
res_a: input array
c1: row index 1
c2: row index 2
"""
res_a = res_a.astype('D')
res_a[r1,:],res_a[r2,:] = scipy.linalg.fblas.zswap(res_a[r1,:],res_a[r2,:])
return (res_a)
def MatrixFloatFix(self,mat,val=1.e-15):
"""
MatrixFloatFix(mat,val=1.e-15)
Clean an array removing any floating point artifacts defined as being smaller than a specified value.
Processes an array inplace
Arguments:
mat: the input 2D array
val [default=1.e-15]: the threshold value (effective zero)
"""
zero_vals = (abs(mat) < val)
for x in range(mat.shape[0]):
for y in range(mat.shape[1]):
if zero_vals[x,y]:
mat[x,y] = 0.0
## if abs(mat[x,y]) != 0.0 and abs(mat[x,y]) < val:
## mat[x,y] = round(mat[x,y])
## if abs(mat[x,y]) < val:
## #mat[x,y] = round(mat[x,y])
## mat[x,y] = 0.0
del zero_vals
def MatrixValueCompare(self,matrix):
"""
MatrixValueCompare(matrix)
Finds the largest/smallest abs(value) > 0.0 in a matrix.
Returns a tuple containing (smallest,largest) values
Arguments:
matrix: the input 2D array
"""
val_B = 0.0
val_S = 1.0e30
for x in range(matrix.shape[0]):
for y in range(matrix.shape[1]):
if abs(matrix[x,y]) > abs(val_B) and abs(matrix[x,y]) != 0.0:
val_B = matrix[x,y]
if abs(matrix[x,y]) < abs(val_S) and abs(matrix[x,y]) != 0.0:
val_S = matrix[x,y]
return (val_S,val_B)
class Stoich(MathArrayFunc):
'''PySCeS stoichiometric analysis class: initialized with a stoichiometric matrix N (input)'''
__stoichdiagmode__ = 0
__version__ = __version__
__TimeFormat = "%H:%M:%S"
USE_QR = False
info_moiety_conserve = False
def __init__(self,input):
"""Initialize class variables"""
self.nmatrix = input
row,col = self.nmatrix.shape
self.nmatrix_col = tuple(scipy.array(range(col)))
self.nmatrix_row = tuple(scipy.array(range(row)))
#Create a machine specific instance
from scipy import MachAr
mach_spec = MachAr()
self.stoichiometric_analysis_fp_zero = mach_spec.eps*2.0e4
self.stoichiometric_analysis_lu_precision = self.stoichiometric_analysis_fp_zero
self.stoichiometric_analysis_gj_precision = self.stoichiometric_analysis_lu_precision*10.0
self.species = None
self.reactions = None
def AnalyseK(self):
"""
AnalyseK()
Evaluate the stoichiometric matrix and calculate the nullspace using LU decomposition and backsubstitution .
Generates the MCA K and Ko arrays and associated row and column vectors
Arguments:
None
"""
print 'Calculating K matrix .',
self.info_flux_conserve = 0 #added 20020416 for conservation detection
if self.__stoichdiagmode__:
print '\nKMATRIX\nGetUpperMatrix: ' + time.strftime(self.__TimeFormat)
p,u,row_vector,column_vector = self.GetUpperMatrix(self.nmatrix)
#p,u,row_vector,column_vector = self.GetUpperMatrixUsingQR(self.nmatrix)
if self.__stoichdiagmode__:
print '\nKMATRIX\nScalePivots: ' + time.strftime(self.__TimeFormat)
unipiv_a = self.ScalePivots(u)
if self.__stoichdiagmode__:
print '\nKMATRIX\nBackSubstitution: ' + time.strftime(self.__TimeFormat)
R_a,row_vector,column_vector = self.BackSubstitution(unipiv_a,row_vector,column_vector)
if self.__stoichdiagmode__:
print '\nKMATRIX\nK_split_R: ' + time.strftime(self.__TimeFormat)
r_ipart,r_fpart,row_vector,column_vector,nullspace,r_fcolumns,self.info_flux_conserve = self.K_split_R(R_a,row_vector,column_vector)
# Don't need anymore ... I think - brett 20051013
## try:
## self.MatrixFloatFix(nullspace,val=self.stoichiometric_analysis_lu_precision)
## except Exception, e:
## if self.__stoichdiagmode__:
## print 'Ignored (no K)',e # brett 20050801
if self.__stoichdiagmode__:
print '\nKMATRIX\n' + time.strftime(self.__TimeFormat)
self.kmatrix = nullspace
self.kmatrix_row = tuple(row_vector)
self.kmatrix_col = tuple(column_vector)
self.kzeromatrix = r_fpart
self.kzeromatrix_row = tuple(r_fcolumns)
self.kzeromatrix_col = tuple(column_vector)
print ' done.'
def AnalyseL(self):
"""
AnalyseL()
Evaluate the stoichiometric matrix and calculate the left nullspace using LU factorization and backsubstitution.
Generates the MCA L, Lo, Nr and Conservation matrix and associated row and column vectors
Arguments:
None
"""
print 'Calculating L matrix .',
a = scipy.transpose(self.nmatrix)
self.info_moiety_conserve = False #added 20020416 for conservation detection
if self.__stoichdiagmode__:
print '\nLMATRIX\nGetUpperMatrix: ' + time.strftime(self.__TimeFormat)
if not self.USE_QR:
p,u,row_vector,column_vector = self.GetUpperMatrix(a)
else:
p,u,row_vector,column_vector = self.GetUpperMatrixUsingQR(a)
if self.__stoichdiagmode__:
print '\nLMATRIX\nScalePivots: ' + time.strftime(self.__TimeFormat)
unipiv_a = self.ScalePivots(u)
if self.__stoichdiagmode__:
print '\nLMATRIX\nBackSubstitution: ' + time.strftime(self.__TimeFormat)
R_a,row_vector,column_vector = self.BackSubstitution(unipiv_a,row_vector,column_vector)
if self.__stoichdiagmode__:
print '\nLMATRIX\nK_split_R: ' + time.strftime(self.__TimeFormat)
r_ipart,consmatrix,cons_row_vector,cons_col_vector,lmatrix,lmatrix_row_vector,\
lmatrix_col_vector,lomatrix,lomatrix_row_vector,lomatrix_col_vector,nrmatrix,Nred_vector,\
Nred_vector_col,self.info_moiety_conserve = self.L_split_R(a,R_a,row_vector,column_vector)
# Don't need anymore ... i think - brett 20051013
## try:
## self.MatrixFloatFix(consmatrix,val=self.stoichiometric_analysis_lu_precision)
## except Exception, e:
## if self.__stoichdiagmode__:
## print 'Ignored (no L)',e # brett 20050801
if self.__stoichdiagmode__:
print '\nLMATRIX\n' + time.strftime(self.__TimeFormat)
self.lmatrix = lmatrix
self.lmatrix_row = tuple(lmatrix_row_vector)
self.lmatrix_col = tuple(lmatrix_col_vector)
self.lzeromatrix = lomatrix
self.lzeromatrix_row = tuple(lomatrix_row_vector)
self.lzeromatrix_col = tuple(lomatrix_col_vector)
self.conservation_matrix = consmatrix
self.conservation_matrix_row = tuple(cons_row_vector)
self.conservation_matrix_col = tuple(cons_col_vector)
self.nrmatrix = nrmatrix
self.nrmatrix_row = tuple(Nred_vector)
self.nrmatrix_col = tuple(scipy.copy(self.nmatrix_col))
print ' done.'
def PivotSort(self,a,row_vector,column_vector):
"""
PivotSort(a,row_vector,column_vector)
This is a sorting routine that accepts a matrix and row/colum vectors
and then sorts them so that: there are no zero rows (by swapping with first
non-zero row) The abs(largest) pivots are moved onto the diagonal to maintain
numerical stability. Row and column swaps are recorded in the tracking vectors.
Arguments:
a: the input array
row_vector: row tracking vector
column_vector: column tracking vector
"""
t = self.commonType(a)
row,col = a.shape
for z in range(0,min(row,col)):
if abs(a[z,z]) < self.stoichiometric_analysis_lu_precision:
maxV = self.stoichiometric_analysis_lu_precision
maxP = (None,None)
zeroP = (None,None)
mList = []
for x in range(z,min(row,col)):
mList.append((x,max(abs(a[x,x:]))))
mVal = 0.0
mRow = None
#print mList
for el in mList:
if el[1] > mVal:
mVal = el[1]
mRow = el[0]
if self.__stoichdiagmode__:
print '\tmVal', mVal, mRow
if mRow != None:
for y in range(z,col):
if abs(a[mRow,y]) > abs(maxV) and abs(a[mRow,y]) > self.stoichiometric_analysis_lu_precision:
maxV = a[mRow,y]
maxP = (mRow,y)
if zeroP != (None,None) and zeroP != (z,z) and mRow != None:
if self.__stoichdiagmode__:
print ' Swapping: ', (z,z), (zeroP[0],zeroP[1]), 1.0
a = self.SwapRowd(a,z,zeroP[0])
a = self.SwapCold(a,z,zeroP[1])
row_vector = self.SwapElem(row_vector,z,zeroP[0])
column_vector = self.SwapElem(column_vector,z,zeroP[1])
elif maxP != (None,None) and maxP != (min(row,col),min(row,col)) and mRow != None:
if self.__stoichdiagmode__:
print ' Swapping: ', (z,z), (maxP[0],maxP[1]), a[z,z], maxV
a = self.SwapRowd(a,z,maxP[0])
a = self.SwapCold(a,z,maxP[1])
row_vector = self.SwapElem(row_vector,z,maxP[0])
column_vector = self.SwapElem(column_vector,z,maxP[1])
return(a,row_vector,column_vector)
def PivotSort_initial(self,a,row_vector,column_vector):
"""
PivotSort_initial(a,row_vector,column_vector)
This is a sorting routine that accepts a matrix and row/colum vectors
and then sorts them so that: the abs(largest) pivots are moved onto the diagonal to maintain
numerical stability i.e. the matrix diagonal is in descending max(abs(value)).
Row and column swaps are recorded in the tracking vectors.
Arguments:
a: the input array
row_vector: row tracking vector
column_vector: column tracking vector
"""
# SAME AS THE ABOVE JUST DOES ALL VALUES NOT ONLY NON_ZERO ONES
# brett 200500802
t = self.commonType(a)
row,col = a.shape
for z in range(0,min(row,col)):
if z == 0 or abs(a[z,z]) < abs(a[z-1,z-1]):
maxV = self.stoichiometric_analysis_lu_precision
maxP = (None,None)
zeroP = (None,None)
mList = []
for x in range(z,min(row,col)):
mList.append((x,max(abs(a[x,x:]))))
mVal = 0.0
mRow = None
for el in mList:
if el[1] > mVal:
mVal = el[1]
mRow = el[0]
if self.__stoichdiagmode__:
print '\tmVal', mVal, mRow
if mRow != None:
for y in range(z,col):
if abs(a[x,y]) > abs(maxV) and abs(a[x,y]) > self.stoichiometric_analysis_lu_precision:
maxV = a[x,y]
maxP = (x,y)
if maxP != (None,None) and maxP != (min(row,col),min(row,col)) and maxP != (z,z):
if self.__stoichdiagmode__:
print ' Swapping: ', (z,z), (maxP[0],maxP[1]), a[z,z], maxV
a = self.SwapRowd(a,z,maxP[0])
a = self.SwapCold(a,z,maxP[1])
row_vector = self.SwapElem(row_vector,z,maxP[0])
column_vector = self.SwapElem(column_vector,z,maxP[1])
return(a,row_vector,column_vector)
def PLUfactorize(self,a_in):
"""
PLUfactorize(a_in)
Performs an LU factorization using LAPACK D/ZGetrf.
Returns LU - combined factorization, IP - rowswap information and info - Getrf error control.
Arguments:
a_in: the matrix to be factorized
"""
print '.',
self.assertRank2(a_in)
t = self.commonType(a_in)
if a_in.dtype.char == 'D':
#print 'Complex matrix ' + a_in.typecode()
a = copy.copy(scipy.transpose(a_in))
elif a_in.dtype.char == 'd':
#print 'Float matrix ' + a_in.typecode()
a = copy.copy(scipy.transpose(a_in))
else:
#print 'Other matrix casting to double, was: ' + a_in.typecode()
a = copy.copy(scipy.transpose(a_in).astype('d'))
Using_FLAPACK = 0
# brett 20041226 - protecting ourselves against flapack
try:
scipy.linalg.clapack.empty_module()
Using_FLAPACK = 1
#print "Using FLAPACK"
except:
pass
#print "Using CLAPACK -- good :-)"
if Using_FLAPACK == 1:
try:
if self.array_kind[t] == 1:
getrf = scipy.linalg.flapack.zgetrf #scipy flapack 20041226
else:
getrf = scipy.linalg.flapack.dgetrf #scipy flapack 20041226
except Exception, e:
print "FLAPACK error", e
else:
try:
if self.array_kind[t] == 1:
getrf = scipy.linalg.clapack.zgetrf #scipy clapack (ATLAS) 20030605
else:
getrf = scipy.linalg.clapack.dgetrf #scipy clapack (ATLAS) 20030605
except Exception, e:
print "CLAPACK error", e
if Using_FLAPACK == 1:
results = getrf(scipy.transpose(a)) # brett 20041226
else:
# This is a $%^&*( ... suddenly with latest cvs f2py getrf only accepts arrays
# not vectors so this song and dance is necessary to fix this 'behaviour'
# as idiotic as it seems, we pad the vector with a zero row (no difference numerically)
# run it through getrf and then strip it afterwards !@#$%^&*() - brett 20050707
if a.shape[0] == 1:
tarr = scipy.zeros((a.shape[0]+1,a.shape[1]),'d')
tarr[0] = a[0]
results = getrf(tarr) #scipy cblas (ATLAS) 20030506 -- this is normal
results = list(results)
results[0] = scipy.array([results[0][0]])
results[1] = scipy.array([results[1][0]],'i')
results[2] = results[2]
results = tuple(results)
else:
results = getrf(a) #scipy cblas (ATLAS) 20030506 -- this is normal
results = list(results)
if results[2] < 0:
print('Argument ', results['info'], ' had an illegal value')
raise LinAlgError
elif results[2] > 0:
pass
if Using_FLAPACK == 1:
result = results[0] # brett 20041226
else:
result = scipy.transpose(results[0]) # -- this is normal
badlist = []
for x in range(result.shape[0]):
for y in range(result.shape[1]):
if abs(result[x,y]) != 0.0 and abs(result[x,y]) < self.stoichiometric_analysis_lu_precision:
result[x,y] = 0.0 # gets rid of the -0.0 irritation
if len(badlist) == 0 and (x,y) == (x,x): # catch 1st float on a pivot
badlist.append(x)
if len(badlist) != 0:
results[2] = badlist[0]-1 # 20040423 if float was on a pivot - refactorize
return(result,results[1],results[2]) #scipy cblas (ATLAS?) 20030506
def SplitLU(self,plu,row,col,t=None):
"""
SplitLU(plu,row,col,t)
PLU takes the combined LU factorization computed by PLUfactorize and extracts the upper matrix.
Returns U.
Arguments:
plu: LU factorization
row: row tracking vector
col: column tracking vector
t [default=None)]: typecode argument (currently not used)
"""
print '.',
if self.__stoichdiagmode__:
print '\n',row,col
for cl in range(0,min(row,col)):
plu[cl+1:,cl] = 0.0
return plu
def GetUpperMatrix(self,a):
"""
GetUpperMatrix(a)
Core analysis algorithm; an input is preconditioned using PivotSort_initial and then cycles of PLUfactorize and
PivotSort are run until the factorization is completed. During this process the matrix is reordered by
column swaps which emulates a full pivoting LU factorization. Returns the pivot matrix P, upper factorization U
as well as the row/col tracking vectors.
Arguments:
a: a stoichiometric matrix
"""
print '.',
t = self.commonType(a)
row,col = a.shape
# this is a test brett 20050802
row_vector = scipy.array((range(row)))
column_vector = scipy.array((range(col)))
a,row_vector,column_vector = self.PivotSort_initial(a,row_vector,column_vector)
#18/09/2000 PLUfactorize included in self.GetUpperMatrix
a,ip,info = self.PLUfactorize(a)
# get U from getrf's PLU
if self.__stoichdiagmode__:
print info
upper_out = self.SplitLU(a,row,col,t)
p_out = scipy.identity(row)
for x in range(0,min(row,col)):
if x + 1 != ip[x]:
p_out = self.SwapRowd(p_out,x,(ip[x]-1))
row_vector = self.SwapElem(row_vector,x,(ip[x]-1))
'''31/08/2000 Added to try to make sure that pivots exist only on the upper left side of the matrix
max(abs(upper_out[x,:])) used instead of abs(max(upper_out[x,:])) otherwise for rows where
all elements < 0 zero is maximum, this way the max of positive values is used'''
upper_out,row_vector,column_vector = self.PivotSort(upper_out,row_vector,column_vector)
'''20/09/2000 This bit sorts out the echelon matrix by running (if needed) cycles of
self.PLUfactorize, self.GetUpperMatrix, and pivsort until only a staircase matrix remains. It uses both the error
generated by self.PLUfactorize(info) and go_flag to control itself'''
if info > 0:
#Here we check to see if the error is not in the last or reduced-last column if it is - exit
# go_flag = 'yes' # 2001/04/26 changed for Python21 and future compatibility
go_flag = 1
for x in range (0,min(row,col)):
if abs(upper_out[x,x]) > self.stoichiometric_analysis_lu_precision:
pos_holder = x
if pos_holder + 1 < info and pos_holder + 1 < row:
if max(abs(upper_out[pos_holder + 1,:])) < self.stoichiometric_analysis_lu_precision:
# go_flag = 'no' # 2001/04/26 changed for Python21 and future compatibility
go_flag = 0
elif pos_holder + 1 == info and pos_holder + 1 == min(row,col):
# go_flag = 'no' # 2001/04/26 changed for Python21 and future compatibility
go_flag = 0
# while info > 0 and go_flag == 'yes': # 2001/04/26 changed for Python21 and future compatibility
while info > 0 and go_flag == 1:
#sort
#upper_out,row_vector,column_vector = pivot_sort2k5(upper_out,row_vector,column_vector)
upper_out,row_vector,column_vector = self.PivotSort(upper_out,row_vector,column_vector)
#PLUfactorize
if self.__stoichdiagmode__:
print ' Echelon: ' + time.strftime(self.__TimeFormat), '(', info,')'
upper_out,ip,info = self.PLUfactorize(upper_out)
# get U from getrf's PLU
if self.__stoichdiagmode__:
print info
upper_out = self.SplitLU(upper_out,row,col,t)
#realign row vector and permutation matrix if necessary
for x in range(0,min(row,col)):
if x + 1 != ip[x]:
p_out = self.SwapRowd(p_out,x,(ip[x]-1))
row_vector = self.SwapElem(row_vector,x,(ip[x]-1))
#test if we can exit -- same criteria as before
for x in range (0,min(row,col)):
if abs(upper_out[x,x]) > self.stoichiometric_analysis_lu_precision:
pos_holder = x
if pos_holder + 1 < info and pos_holder + 1 < row:
if max(abs(upper_out[pos_holder + 1,:])) < self.stoichiometric_analysis_lu_precision:
# go_flag = 'no' # 2001/04/26 changed for Python21 and future compatibility
go_flag = 0
elif pos_holder + 1 == info and pos_holder + 1 == min(row,col):
# go_flag = 'no' # 2001/04/26 changed for Python21 and future compatibility
go_flag = 0
'''This bit will get rid of any zero rows so that we will hopefully only have to work with a
reduced matrix after this, for completeness the row_vector will also be sliced'''
for x in range (0,min(row,col)):
if abs(upper_out[x,x]) > self.stoichiometric_analysis_lu_precision:
pos_holder = x
upper_out_r = scipy.zeros((pos_holder + 1,col)).astype(t)
upper_out_r = upper_out[:pos_holder + 1,:]
row_vector_r = row_vector[:pos_holder + 1]
return(p_out,upper_out_r,row_vector_r,column_vector)
def GetUpperMatrixUsingQR(self,a):
"""
GetUpperMatrix(a)
Core analysis algorithm; an input is preconditioned using PivotSort_initial and then cycles of PLUfactorize and
PivotSort are run until the factorization is completed. During this process the matrix is reordered by
column swaps which emulates a full pivoting LU factorization. Returns the pivot matrix P, upper factorization U
as well as the row/col tracking vectors.
Arguments:
a: a stoichiometric matrix
"""
print '.',
t = self.commonType(a)
row,col = a.shape
# this is a test brett 20050802
row_vector = scipy.array((range(row)))
column_vector = scipy.array((range(col)))
## a,row_vector,column_vector = self.PivotSort(a,row_vector,column_vector)
a,row_vector,column_vector = self.PivotSort_initial(a,row_vector,column_vector)
Q,upper_out = scipy.linalg.qr(a)
self.MatrixFloatFix(upper_out,val=self.stoichiometric_analysis_lu_precision*10.0)
'''This bit will get rid of any zero rows so that we will hopefully only have to work with a
reduced matrix after this, for completeness the row_vector will also be sliced'''
for x in range (0,min(row,col)):
if abs(upper_out[x,x]) > self.stoichiometric_analysis_lu_precision:
pos_holder = x
upper_out_r = scipy.zeros((pos_holder + 1,col)).astype(t)
p_out = scipy.diag(upper_out.shape[1]*[1.0])
upper_out_r = upper_out[:pos_holder + 1,:]
row_vector_r = row_vector[:pos_holder + 1]
return(p_out,upper_out_r,row_vector_r,column_vector)
def ScalePivots(self,a_one):
"""
ScalePivots(a_one)
Given an upper triangular matrix U, this method scales the diagonal (pivot values) to one.
Arguments:
a_one: an upper triangular matrix U
"""
print '.',
t = self.commonType(a_one)
row,col = a_one.shape
'''13/09/2000 We now assume that the matrix has the correct shape ie. the pivots are in a
perfect staircase'''
for x in range (0,row):
if abs(a_one[x,x]) == 0.0: # ths is to prevent NAN's when FixFloat zeros a supersmall pivot
pass
elif abs(a_one[x,x]) != 1.0:
a_one[x,:] = a_one[x,:]/abs(a_one[x,x])
if a_one[x,x] < 0.0:
a_one[x,:] = -(a_one[x,:])
return(a_one)
def BackSubstitution(self,res_a,row_vector,column_vector):
"""
BackSubstitution(res_a,row_vector,column_vector)
Jordan reduction of a scaled upper triangular matrix. The returned array is now in the form [I R] and can
be used for nullspace determination. Modified row and column tracking vetors are also returned.
Arguments:
res_a: unitary pivot upper triangular matrix
row_vector: row tracking vector
column_vector: column tracking vector
"""
print '.',
t = self.commonType(res_a)
row,col = res_a.shape
'''13/09/2000 removed the copy thing, and eliminated the divnumb variable. Because we
are now working with a row reduced matrix upward elimination only works on the pivot cols'''
# old back substitution circa 2000! brett - 20050801
## bigF = 0
## if max(res_a.shape) > 500:
## bigF = 1
## for y in range (min(row,col)-1,-1,-1):
## if bigF and self.__stoichdiagmode__:
## print '.',
## for x in range (min(row,col)-2,-1,-1):
## z = x + 1
## while z < min(row,col):
## if abs(res_a[x,y]) > self.stoichiometric_analysis_gj_precision and abs(res_a[z,y]) > self.stoichiometric_analysis_gj_precision:
## res_a[x,:] = res_a[x,:] - res_a[z,:]*(res_a[x,y]/res_a[z,y])
## z = z + 1
## continue
## else:
## z = z + 1
# new back substitution
# right looking algorithm that uses array slicing speeds up as it moves down the pivots
# this algorithm is at least 3X as fast as the old one for a large system - brett 20050805
bigF = 0
if max(res_a.shape) > 500:
bigF = 1
for x in range(min(row,col)):
if bigF and self.__stoichdiagmode__:
print '.',
for y in range(x+1,min(row,col)):
if abs(res_a[x,y]) > self.stoichiometric_analysis_lu_precision:
res_a[x,y:] = res_a[x,y:]-(res_a[x,y]*res_a[y,y:])
# $%^&* wtf is this for??? brett - 20050801 (besides cleaning fp wierdness ... nothing ;-)
self.MatrixFloatFix(res_a,val=self.stoichiometric_analysis_gj_precision)
res_a,row_vector,column_vector = self.PivotSort(res_a,row_vector,column_vector)
return (res_a,row_vector,column_vector)
def K_split_R(self,R_a,row_vector,column_vector):
"""
K_split_R(R_a,row_vector,column_vector)
Using the R factorized form of the stoichiometric matrix we now form the K and Ko matrices. Returns
the r_ipart,Komatrix,Krow,Kcolumn,Kmatrix,Korow,info
Arguments:
R_a: the Gauss-Jordan reduced stoichiometric matrix
row_vector: row tracking vector
column_vector: column tracking vector
"""
print '.',
t = self.commonType(R_a)
row,col = R_a.shape
'''14/09/2000 Seeing as we now should have a perfect staircase in a reduced matrix we do
not have to search for the last pivot it will exist at min(row,col)-1 so the pos_holder
finding code has been removed (actually moved into self.GetUpperMatrix)'''
pos_holder = min(row,col) - 1
'''This bit extracts the identity part from R (future note this could be replaced by an I matrix formed by min(row,col))'''
r_ipart = scipy.zeros((pos_holder + 1,pos_holder + 1)).astype(t)
r_ipart = R_a[:pos_holder+1,:pos_holder+1]
row_i,col_i = r_ipart.shape
'''If there are free variables, then this bit will extract them and form the row/col vectors'''
empty_rf = 0
K_switch = 0 #added 20020416 class attribute for conservation detection 0 = none, 1 = exists
if col - col_i > self.stoichiometric_analysis_fp_zero:
r_fpart = scipy.zeros((pos_holder + 1,col - (pos_holder + 1))).astype(t)
r_fpart = R_a[:pos_holder+1,pos_holder+1:]
row_vector_dependent = column_vector[:pos_holder + 1]
row_vector_independent = column_vector[pos_holder + 1:]
row_vector = scipy.concatenate((row_vector_independent,row_vector_dependent),1)
column_vector = column_vector[pos_holder + 1:]
K_switch = 1
else:
print('no flux conservation')
r_fpart = 'no flux conservation'
empty_rf = 1
K_switch = 0
# if r_fpart == 'F is empty, R = I': # 2001/04/26 changed for Python21 compatibility
if empty_rf == 1:
return(r_ipart,r_ipart,column_vector,column_vector,r_ipart,column_vector,K_switch)
else:
row,col = r_fpart.shape
id = scipy.identity(col).astype(t)
nullspace = scipy.concatenate((id,-r_fpart),0).astype(t)
#brett 05/11/2002 changed r_fpart to -r_fpart
return(r_ipart,-r_fpart,row_vector,column_vector,nullspace,row_vector_dependent,K_switch)
def L_split_R(self,Nfull,R_a,row_vector,column_vector):
"""
L_split_R(Nfull,R_a,row_vector,column_vector)
Takes the Gauss-Jordan factorized N^T and extract the L, Lo, conservation (I -Lo) and reduced stoichiometric matrices. Returns: lmatrix_col_vector, lomatrix, lomatrix_row, lomatrix_co, nrmatrix, Nred_vector_row, Nred_vector_col, info
Arguments:
Nfull: the original stoichiometric matrix N
R_a: gauss-jordan factorized form of N^T
row_vector: row tracking vector
column_vector: column tracking vector
"""
print '.',
#print '\nR_a'
#print `R_a`
t = self.commonType(R_a)
row,col = R_a.shape
'''14/09/2000 Seeing as we now should have a perfect staircase in a reduced matrix we do
not have to search for the last pivot it will exist at min(row,col)-1 so the pos_holder
finding code has been removed (actually moved into self.GetUpperMatrix)'''
pos_holder = min(row,col) - 1
'''Here we extract the identity matrix from R (future note this could be replaced by an I matrix formed by min(row,col))'''
r_ipart = scipy.zeros((pos_holder + 1,pos_holder + 1)).astype(t)
r_ipart = R_a[:pos_holder+1,:pos_holder+1]
row_i,col_i = r_ipart.shape
# exit1 = 'no' # 2001/04/26 changed for Python21 and future compatibility
exit1 = 0
L_switch = False #added 20020416 class attribute for conservation detection 0 = none, 1 = exists
'''If there are free variable then they are extracted and packaged, row/col vectors are formed'''
if col - col_i > self.stoichiometric_analysis_fp_zero:
r_fpart = scipy.zeros((pos_holder + 1,col - (pos_holder + 1))).astype(t)
r_fpart = R_a[:pos_holder+1,pos_holder+1:]
col_vector_dependent = column_vector[pos_holder + 1:]
col_vector_independent = column_vector[:pos_holder + 1]
cons_col_vector = scipy.concatenate((col_vector_independent,col_vector_dependent),1)
cons_row_vector = col_vector_dependent
lmatrix_row_vector = scipy.concatenate((col_vector_independent,col_vector_dependent),1)
lmatrix_col_vector = col_vector_independent
lomatrix_row_vector = col_vector_dependent
lomatrix_col_vector = col_vector_independent
Nred_vector = col_vector_independent
L_switch = True
else:
lomatrix_row_vector = column_vector[:pos_holder + 1]
lomatrix_col_vector = column_vector[:pos_holder + 1]
Nred_vector = column_vector[:pos_holder + 1]
#print('F is empty, R = I, no conservation matrix, L is Lo')
r_fpart = r_ipart
# exit1 = 'yes' # 2001/04/26 changed for Python21 and future compatibility
exit1 = 1
L_switch = False
# if exit1 == 'yes': # 2001/04/26 changed for Python21 and future compatibility
if exit1 == 1:
r_fpart = scipy.transpose(r_fpart)
lmatrix = r_fpart
#02/10/2000 removed so that the thing returns the Lo matrix as L
#r_fpart = 'no conservation Lo = L = I'
# my factorization routines now swap things around for numeric stability, so that Nr might be a row/column swapped
# this simply synchronizes Nr with its labels - brett 20050805
Nfull = scipy.transpose(Nfull)
row,col = Nfull.shape
Nred = scipy.zeros((row_i,col)).astype(t)
for x in range (0,row_i):
Nred[x,:] = Nfull[Nred_vector[x],:]
# return the right stuff - brett 20050805
return(r_ipart,'no conservation','no conservation','no conservation',lmatrix,lomatrix_row_vector,lomatrix_col_vector,lmatrix,lomatrix_row_vector,lomatrix_col_vector,Nred,Nred_vector,row_vector,L_switch)
#return(r_ipart,'no conservation','no conservation','no conservation',lmatrix,lomatrix_row_vector,lomatrix_col_vector,lmatrix,lomatrix_row_vector,lomatrix_col_vector,scipy.transpose(Nfull),Nred_vector,row_vector,L_switch)
else:
r_fpart = scipy.transpose(r_fpart)
row,col = r_fpart.shape
id = scipy.identity(row).astype(t)
consmatrix = scipy.concatenate((-r_fpart,id),1).astype(t)
id = scipy.identity(col).astype(t)
lmatrix = scipy.concatenate((id,r_fpart),0).astype(t)
'''This bit creates Nr. The transpose is only necessary if Nfull is already transposed in the input function'''
Nfull = scipy.transpose(Nfull)
row,col = Nfull.shape
Nred = scipy.zeros((row_i,col)).astype(t)
for x in range (0,row_i):
Nred[x,:] = Nfull[Nred_vector[x],:]
return(r_ipart,consmatrix,cons_row_vector,cons_col_vector,lmatrix,lmatrix_row_vector,lmatrix_col_vector,r_fpart,lomatrix_row_vector,lomatrix_col_vector,Nred,Nred_vector,row_vector,L_switch)
def SVD_Rank_Check(self,matrix=None,factor=1.0e4,resultback=0):
"""
SVD_Rank_Check(matrix=None,factor=1.0e4,resultback=0)
Calculates the dimensions of L/L0/K/K) by way of SVD and compares them to the Guass-Jordan results. Please note that for LARGE ill conditioned matrices the SVD can become numerically unstable when used for nullspace determinations
Arguments:
matrix [default=None]: the stoichiometric matrix default is self.Nmatrix
factor [default=1.0e4]: factor used to calculate the 'zero pivot' mask = mach_eps*factor
resultback [default=0]: return the SVD results, U, S, vh
"""
if matrix == None:
matrix = self.nmatrix
nrow,ncol = matrix.shape
TrMat = 0
if ncol > nrow:
# 'INFO: SVD is more accurate for tall matrices using (N)T\n'
matrix = scipy.transpose(matrix)
TrMat = 1
u,s,vh = scipy.linalg.svd(matrix)
maskF = scipy.machar.machar_double.eps*factor
if TrMat == 0:
print 'SVD zero mask:', maskF
print 'LU factorization effective zero:', self.stoichiometric_analysis_lu_precision
rank = len([el for el in (abs(s) > maskF) if el > 0.0])
## null_mask = (abs(s) > maskF)
## vhnull = scipy.compress(null_mask,vh,axis=0)
## unull = scipy.compress(null_mask,u,axis=0)
## assert vhnull.shape[0] == unull.shape[0], 'This should be true or I\'m very confused'
## rank = vhnull.shape[0]
print '\nNmatrix has ' + `nrow` + ' rows and ' + `ncol` + ' columns'
print '\nSVD \"considers\" the rank to be: ' + `rank`
print 'LU (Kmatrix) considers the rank to be: ' + `self.kzeromatrix.shape[0]`
print 'LU (Lmatrix) considers the rank to be: ' + `self.lzeromatrix.shape[1]`
print '\nComparing lzeromatrix dimensions'
print 'SVD lzeromatrix.shape =', (u.shape[0]-rank,rank)
print 'LU lzeromatrix.shape =', self.lzeromatrix.shape
print '\nComparing lmatrix dimensions'
print 'SVD lmatrix.shape =', (u.shape[0],rank)
print 'LU lmatrix.shape =', self.lmatrix.shape
print '\nComparing kzeromatrix dimensions'
print 'SVD kzeromatrix.shape =', (rank,vh.shape[0]-rank)
print 'LU kzeromatrix.shape =', self.kzeromatrix.shape
print '\nComparing kmatrix dimensions'
print 'SVD kmatrix.shape =', (vh.shape[0],vh.shape[0]-rank)
print 'LU kmatrix.shape =', self.kmatrix.shape
else:
print 'SVD zero mask:', maskF
print 'LU factorization effective zero:', self.stoichiometric_analysis_lu_precision
rank = len([el for el in (abs(s) > maskF) if el > 0.0])
## null_mask = (abs(s) > maskF)
## vhnull = scipy.compress(null_mask,vh,axis=0)
## unull = scipy.compress(null_mask,u,axis=0)
## assert vhnull.shape[0] == unull.shape[0], 'This should be true or I\'m very confused'
## rank = vhnull.shape[0]
print '\nNmatrix has ' + `nrow` + ' rows and ' + `ncol` + ' columns'
print '\nSVD considers the rank to be: ' + `rank`
print 'LU (Kmatrix) considers the rank to be: ' + `self.kzeromatrix.shape[0]`
print 'LU (Lmatrix) considers the rank to be: ' + `self.lzeromatrix.shape[1]`
print '\nComparing lzeromatrix dimensions'
print 'SVD lzeromatrix.shape =', (vh.shape[0]-rank,rank)
print 'LU lzeromatrix.shape =', self.lzeromatrix.shape
print '\nComparing lmatrix dimensions'
print 'SVD lmatrix.shape =', (vh.shape[0],rank)
print 'LU lmatrix.shape =', self.lmatrix.shape
print '\nComparing kzeromatrix dimensions'
print 'SVD kzeromatrix.shape =', (rank,u.shape[0]-rank)
print 'LU kzeromatrix.shape =', self.kzeromatrix.shape
print '\nComparing kmatrix dimensions'
print 'SVD kmatrix.shape =', (u.shape[0],u.shape[0]-rank)
print 'LU kmatrix.shape =', self.kmatrix.shape
print '\nMatrix value check\n******************'
print '\nKzeromatrix:'
sm,bg = self.MatrixValueCompare(self.kzeromatrix)
print 'Smallest value: ', abs(sm)
print 'Largest value: ', abs(bg)
print 'Ratio abs(bg/sm):', abs(bg/sm)
print '\nLzeromatrix:'
sm,bg = self.MatrixValueCompare(self.lzeromatrix)
print 'Smallest value: ', abs(sm)
print 'Largest value: ', abs(bg)
print 'Ratio abs(bg/sm):', abs(bg/sm)
print '\nPlease note: I\'ve found that for larger models the rank calculated by SVD in this test can become unstable and give incorrect results.'
print ' '
if resultback:
return u,s,vh
if __psyco_active__:
psyco.bind(MathArrayFunc)
psyco.bind(Stoich)
|