1""" 2PySCeS - Python Simulator for Cellular Systems (http://pysces.sourceforge.net) 3 4Copyright (C) 2004-2020 B.G. Olivier, J.M. Rohwer, J.-H.S Hofmeyr all rights reserved, 5 6Brett G. Olivier (bgoli@users.sourceforge.net) 7Triple-J Group for Molecular Cell Physiology 8Stellenbosch University, South Africa. 9 10Permission to use, modify, and distribute this software is given under the 11terms of the PySceS (BSD style) license. See LICENSE.txt that came with 12this distribution for specifics. 13 14NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK. 15Brett G. Olivier 16""" 17from __future__ import division, print_function 18from __future__ import absolute_import 19from __future__ import unicode_literals 20 21from pysces.version import __version__ 22import time,os 23import scipy 24import scipy.linalg 25 26# if int(scipy.__version__.split('.')[1]) < 12: 27# thanks scipy 28if int(scipy.__version__.split('.')[0]) < 1 and int(scipy.__version__.split('.')[1]) < 12: 29 myfblas = scipy.linalg.fblas 30 myflapack = scipy.linalg.flapack 31else: 32 myfblas = scipy.linalg.blas 33 myflapack = scipy.linalg.lapack 34 35import copy 36 37__doc__ = """ 38 PyscesStoich 39 ------------ 40 PySCeS stoichiometric analysis classes. 41 42 """ 43 44## print 'Stoichiometry ver ' + __version__ + ' runtime: '+ time.strftime("%H:%M:%S") 45 46__psyco_active__ = 0 47## try: 48 ## import psyco 49 ## psyco.profile() 50 ## __psyco_active__ = 1 51 ## print 'PySCeS Stoichiometry Module is now PsycoActive!' 52## except: 53 ## __psyco_active__ = 0 54 55##stoich_zero_valM = scipy.machar.MachAr().eps*2.0e4 # safer --> about 4e-11 (unofficially - a minpivot size) 56##stoich_zero_valM = scipy.machar.MachAr().eps*2.0e4 # safe --> about 4e-12 (unofficially - a minpivot size) 57##print '\tStoichiometric precision = ', stoich_zero_valM 58 59class StructMatrix: 60 """ 61 This class is specifically designed to store structural matrix information 62 give it an array and row/col index permutations it can generate its own 63 row/col labels given the label src. 64 """ 65 66 array = None 67 ridx = None 68 cidx = None 69 row = None 70 col = None 71 shape = None 72 73 def __init__(self, array, ridx, cidx, row=None, col=None): 74 """ 75 Instantiate with array and matching row/col index arrays, optional label arrays 76 """ 77 self.array = array 78 self.ridx = ridx 79 self.cidx = cidx 80 self.row = row 81 self.col = col 82 self.shape = array.shape 83 84 def __call__(self): 85 return self.array 86 87 def getRowsByIdx(self, *args): 88 """Return the rows referenced by index (1,3,5)""" 89 return self.array.take(args, axis=0) 90 91 def getColsByIdx(self, *args): 92 """Return the columns referenced by index (1,3,5)""" 93 return self.array.take(args, axis=1) 94 95 def setRow(self, src): 96 """ 97 Assuming that the row index array is a permutation (full/subset) 98 of a source label array by supplying that source to setRow it 99 maps the row labels to ridx and creates self.row (row label list) 100 """ 101 self.row = [src[r] for r in self.ridx] 102 103 def setCol(self, src): 104 """ 105 Assuming that the col index array is a permutation (full/subset) 106 of a source label array by supplying that src to setCol 107 maps the row labels to cidx and creates self.col (col label list) 108 """ 109 110 self.col = [src[c] for c in self.cidx] 111 112 def getRowsByName(self, *args): 113 """Return the rows referenced by label ('s','x','d')""" 114 assert self.row != None, "\nI need row labels" 115 try: 116 return self.array.take([self.row.index(l) for l in args], axis=0) 117 except Exception as ex: 118 print(ex) 119 print("\nValid row labels are: %s" % self.row) 120 return None 121 122 def getColsByName(self, *args): 123 """Return the columns referenced by label ('s','x','d')""" 124 assert self.col != None, "\nI need column labels" 125 try: 126 return self.array.take([self.col.index(l) for l in args], axis=1) 127 except Exception as ex: 128 print(ex) 129 print("Valid column labels are: %s" % self.col) 130 return None 131 132 def getLabels(self, axis='all'): 133 """Return the matrix labels ([rows],[cols]) where axis='row'/'col'/'all'""" 134 if axis == 'row': return self.row 135 elif axis == 'col': return self.col 136 else: return self.row, self.col 137 138 def getIndexes(self, axis='all'): 139 """Return the matrix indexes ([rows],[cols]) where axis='row'/'col'/'all'""" 140 if axis == 'row': return self.ridx 141 elif axis == 'col': return self.cidx 142 else: return self.ridx, self.cidx 143 144 def getByIdx(self, row, col): 145 assert row in self.ridx, '\n%s is an invalid index' % row 146 assert col in self.cidx, '\n%s is an invalid index' % col 147 return self.array[row, col] 148 149 def getByName(self, row, col): 150 assert row in self.row, '\n%s is an invalid name' % row 151 assert col in self.col, '\n%s is an invalid name' % col 152 return self.array[self.row.index(row), self.col.index(col)] 153 154 def setByIdx(self, row, col, val): 155 assert row in self.ridx, '\n%s is an invalid index' % row 156 assert col in self.cidx, '\n%s is an invalid index' % col 157 self.array[row, col] = val 158 159 def setByName(self, row, col, val): 160 assert row in self.row, '\n%s is an invalid name' % row 161 assert col in self.col, '\n%s is an invalid name' % col 162 self.array[self.row.index(row), self.col.index(col)] = val 163 164class MathArrayFunc(object): 165 """A class of basic array functions some LAPACK based""" 166 167 __doc__ = '''PySCeS array functions - used by Stoich''' 168 169 array_kind = {'i':0, 'l': 0, 'f': 0, 'd': 0, 'F': 1, 'D': 1} 170 array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1} 171 array_type = [['f', 'd'], ['F', 'D']] 172 LinAlgError = 'LinearAlgebraError' 173 174 def commonType(self,*arrays): 175 """ 176 commonType(\*arrays) 177 178 Numeric detect and set array precision (will be replaced with new scipy.core compatible code when ready) 179 180 Arguments: 181 182 \*arrays: input arrays 183 184 """ 185 kind = 0 186 # precision = 0 187 # force higher precision in lite version 188 precision = 1 189 for a in arrays: 190 t = a.dtype.char 191 kind = max(kind, self.array_kind[t]) 192 precision = max(precision, self.array_precision[t]) 193 194 return self.array_type[kind][precision] 195 196 def castCopyAndTranspose(self,type, *arrays): 197 """ 198 castCopyAndTranspose(type, \*arrays) 199 200 Cast numeric arrays to required type and transpose 201 202 Arguments: 203 204 type: the required type to cast to 205 \*arrays: the arrays to be processed 206 207 """ 208 cast_arrays = () 209 for a in arrays: 210 if a.typecode() == type: 211 cast_arrays = cast_arrays + (copy.copy(scipy.transpose(a)),) 212 else: 213 cast_arrays = cast_arrays + (copy.copy(scipy.transpose(a).astype(type)),) 214 if len(cast_arrays) == 1: 215 return cast_arrays[0] 216 else: 217 return cast_arrays 218 219 def assertRank2(self,*arrays): 220 """ 221 assertRank2(\*arrays) 222 223 Check that we are using a 2D array 224 225 Arguments: 226 227 \*arrays: input array(s) 228 229 """ 230 for a in arrays: 231 if len(a.shape) != 2: 232 raise LinAlgError('Array must be two-dimensional') 233 234 235 def SwapCol(self,res_a,r1,r2): 236 """ 237 SwapCol(res_a,r1,r2) 238 239 Swap two columns using BLAS swap, arrays can be (or are upcast to) type double (d) or double complex (D). 240 Returns the colswapped array 241 242 Arguments: 243 244 res_a: the input array 245 r1: the first column to be swapped 246 r2: the second column to be swapped 247 248 """ 249 if self.array_kind[self.commonType(res_a)] == 1: # brett 20041226 added complex support to swap 250 res_a = res_a.astype('D') 251 return self.SwapColz(res_a,r1,r2) 252 else: 253 res_a = res_a.astype('d') 254 return self.SwapCold(res_a,r1,r2) 255 256 257 def SwapRow(self,res_a,r1,r2): 258 """ 259 SwapRow(res_a,r1,r2) 260 261 Swaps two rows using BLAS swap, arrays can be (or are upcast to) type double (d) or double complex (D). 262 Returns the rowswapped array. 263 264 Arguments: 265 266 res_a: the input array 267 r1: the first row index to be swapped 268 r2: the second row index to be swapped 269 270 """ 271 if self.array_kind[self.commonType(res_a)] == 1: # brett 20041226 added complex support to swap 272 res_a = res_a.astype('D') 273 return self.SwapRowz(res_a,r1,r2) 274 else: 275 res_a = res_a.astype('d') 276 return self.SwapRowd(res_a,r1,r2) 277 278 def SwapElem(self,res_a,r1,r2): 279 """ 280 SwapElem(res_a,r1,r2) 281 282 Swaps two elements in a 1D vector 283 284 Arguments: 285 286 res_a: the input vector 287 r1: index 1 288 r2: index 2 289 290 """ 291 res_a[r1],res_a[r2] = res_a[r2],res_a[r1] 292 return (res_a) 293 294 def SwapCold(self,res_a,c1,c2): 295 """ 296 SwapCold(res_a,c1,c2) 297 298 Swaps two double (d) columns in an array using BLAS DSWAP. Returns the colswapped array. 299 300 Arguments: 301 302 res_a: input array 303 c1: column index 1 304 c2: column index 2 305 306 """ 307 res_a = res_a.astype('d') 308 res_a[:,c1],res_a[:,c2] = myfblas.dswap(res_a[:,c1],res_a[:,c2]) 309 return (res_a) 310 311 def SwapRowd(self,res_a,r1,r2): 312 """ 313 SwapRowd(res_a,c1,c2) 314 315 Swaps two double (d) rows in an array using BLAS DSWAP. Returns the rowswapped array. 316 317 Arguments: 318 319 res_a: input array 320 c1: row index 1 321 c2: row index 2 322 323 """ 324 res_a = res_a.astype('d') 325 res_a[r1,:],res_a[r2,:] = myfblas.dswap(res_a[r1,:],res_a[r2,:]) 326 return (res_a) 327 328 def SwapColz(self,res_a,c1,c2): 329 """ 330 SwapColz(res_a,c1,c2) 331 332 Swaps two double complex (D) columns in an array using BLAS ZSWAP. Returns the colswapped array. 333 334 Arguments: 335 336 res_a: input array 337 c1: column index 1 338 c2: column index 2 339 340 """ 341 res_a = res_a.astype('D') 342 res_a[:,c1],res_a[:,c2] = myfblas.zswap(res_a[:,c1],res_a[:,c2]) 343 return (res_a) 344 345 def SwapRowz(self,res_a,r1,r2): 346 """ 347 SwapRowz(res_a,c1,c2) 348 349 Swaps two double complex (D) rows in an array using BLAS ZSWAP. Returns the rowswapped array. 350 351 Arguments: 352 353 res_a: input array 354 c1: row index 1 355 c2: row index 2 356 357 """ 358 res_a = res_a.astype('D') 359 res_a[r1,:],res_a[r2,:] = myfblas.zswap(res_a[r1,:],res_a[r2,:]) 360 return (res_a) 361 362 def MatrixFloatFix(self,mat,val=1.e-15): 363 """ 364 MatrixFloatFix(mat,val=1.e-15) 365 366 Clean an array removing any floating point artifacts defined as being smaller than a specified value. 367 Processes an array inplace 368 369 Arguments: 370 371 mat: the input 2D array 372 val [default=1.e-15]: the threshold value (effective zero) 373 374 """ 375 zero_vals = (abs(mat) < val) 376 377 for x in range(mat.shape[0]): 378 for y in range(mat.shape[1]): 379 if zero_vals[x,y]: 380 mat[x,y] = 0.0 381 ## if abs(mat[x,y]) != 0.0 and abs(mat[x,y]) < val: 382 ## mat[x,y] = round(mat[x,y]) 383 ## if abs(mat[x,y]) < val: 384 ## #mat[x,y] = round(mat[x,y]) 385 ## mat[x,y] = 0.0 386 del zero_vals 387 388 def MatrixValueCompare(self,matrix): 389 """ 390 MatrixValueCompare(matrix) 391 392 Finds the largest/smallest abs(value) > 0.0 in a matrix. 393 Returns a tuple containing (smallest,largest) values 394 395 Arguments: 396 397 matrix: the input 2D array 398 399 """ 400 val_B = 0.0 401 val_S = 1.0e30 402 for x in range(matrix.shape[0]): 403 for y in range(matrix.shape[1]): 404 if abs(matrix[x,y]) > abs(val_B) and abs(matrix[x,y]) != 0.0: 405 val_B = matrix[x,y] 406 if abs(matrix[x,y]) < abs(val_S) and abs(matrix[x,y]) != 0.0: 407 val_S = matrix[x,y] 408 return (val_S,val_B) 409 410 411class Stoich(MathArrayFunc): 412 '''PySCeS stoichiometric analysis class: initialized with a stoichiometric matrix N (input)''' 413 __stoichdiagmode__ = 0 414 __version__ = __version__ 415 __TimeFormat = "%H:%M:%S" 416 USE_QR = False 417 info_moiety_conserve = False 418 419 def __init__(self,input): 420 """Initialize class variables""" 421 422 self.nmatrix = input 423 row,col = self.nmatrix.shape 424 self.nmatrix_col = tuple(scipy.array(list(range(col)))) 425 self.nmatrix_row = tuple(scipy.array(list(range(row)))) 426 427 #Create a machine specific instance 428 from scipy import MachAr 429 mach_spec = MachAr() 430 431 self.stoichiometric_analysis_fp_zero = mach_spec.eps*2.0e4 432 self.stoichiometric_analysis_lu_precision = self.stoichiometric_analysis_fp_zero 433 self.stoichiometric_analysis_gj_precision = self.stoichiometric_analysis_lu_precision*10.0 434 435 self.species = None 436 self.reactions = None 437 438 def AnalyseK(self): 439 """ 440 AnalyseK() 441 442 Evaluate the stoichiometric matrix and calculate the nullspace using LU decomposition and backsubstitution . 443 Generates the MCA K and Ko arrays and associated row and column vectors 444 445 Arguments: 446 None 447 448 """ 449 print('Calculating K matrix .', end=' ') 450 451 self.info_flux_conserve = 0 #added 20020416 for conservation detection 452 453 if self.__stoichdiagmode__: 454 print('\nKMATRIX\nGetUpperMatrix: ' + time.strftime(self.__TimeFormat)) 455 456 p,u,row_vector,column_vector = self.GetUpperMatrix(self.nmatrix) 457 #p,u,row_vector,column_vector = self.GetUpperMatrixUsingQR(self.nmatrix) 458 if self.__stoichdiagmode__: 459 print('\nKMATRIX\nScalePivots: ' + time.strftime(self.__TimeFormat)) 460 461 unipiv_a = self.ScalePivots(u) 462 if self.__stoichdiagmode__: 463 print('\nKMATRIX\nBackSubstitution: ' + time.strftime(self.__TimeFormat)) 464 465 R_a,row_vector,column_vector = self.BackSubstitution(unipiv_a,row_vector,column_vector) 466 if self.__stoichdiagmode__: 467 print('\nKMATRIX\nK_split_R: ' + time.strftime(self.__TimeFormat)) 468 469 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) 470 # Don't need anymore ... I think - brett 20051013 471## try: 472## self.MatrixFloatFix(nullspace,val=self.stoichiometric_analysis_lu_precision) 473## except Exception, e: 474## if self.__stoichdiagmode__: 475## print 'Ignored (no K)',e # brett 20050801 476 477 478 if self.__stoichdiagmode__: 479 print('\nKMATRIX\n' + time.strftime(self.__TimeFormat)) 480 481 self.kmatrix = nullspace 482 self.kmatrix_row = tuple(row_vector) 483 self.kmatrix_col = tuple(column_vector) 484 485 self.kzeromatrix = r_fpart 486 self.kzeromatrix_row = tuple(r_fcolumns) 487 self.kzeromatrix_col = tuple(column_vector) 488 print(' done.') 489 490 def AnalyseL(self): 491 """ 492 AnalyseL() 493 494 Evaluate the stoichiometric matrix and calculate the left nullspace using LU factorization and backsubstitution. 495 Generates the MCA L, Lo, Nr and Conservation matrix and associated row and column vectors 496 497 Arguments: 498 None 499 500 """ 501 print('Calculating L matrix .', end=' ') 502 503 a = scipy.transpose(self.nmatrix) 504 self.info_moiety_conserve = False #added 20020416 for conservation detection 505 506 if self.__stoichdiagmode__: 507 print('\nLMATRIX\nGetUpperMatrix: ' + time.strftime(self.__TimeFormat)) 508 509 510 if not self.USE_QR: 511 p,u,row_vector,column_vector = self.GetUpperMatrix(a) 512 else: 513 p,u,row_vector,column_vector = self.GetUpperMatrixUsingQR(a) 514 if self.__stoichdiagmode__: 515 print('\nLMATRIX\nScalePivots: ' + time.strftime(self.__TimeFormat)) 516 517 unipiv_a = self.ScalePivots(u) 518 if self.__stoichdiagmode__: 519 print('\nLMATRIX\nBackSubstitution: ' + time.strftime(self.__TimeFormat)) 520 521 R_a,row_vector,column_vector = self.BackSubstitution(unipiv_a,row_vector,column_vector) 522 if self.__stoichdiagmode__: 523 print('\nLMATRIX\nK_split_R: ' + time.strftime(self.__TimeFormat)) 524 525 r_ipart,consmatrix,cons_row_vector,cons_col_vector,lmatrix,lmatrix_row_vector,\ 526 lmatrix_col_vector,lomatrix,lomatrix_row_vector,lomatrix_col_vector,nrmatrix,Nred_vector,\ 527 Nred_vector_col,self.info_moiety_conserve = self.L_split_R(a,R_a,row_vector,column_vector) 528 # Don't need anymore ... i think - brett 20051013 529 ## try: 530 ## self.MatrixFloatFix(consmatrix,val=self.stoichiometric_analysis_lu_precision) 531 ## except Exception, e: 532 ## if self.__stoichdiagmode__: 533 ## print 'Ignored (no L)',e # brett 20050801 534 535 if self.__stoichdiagmode__: 536 print('\nLMATRIX\n' + time.strftime(self.__TimeFormat)) 537 538 self.lmatrix = lmatrix 539 self.lmatrix_row = tuple(lmatrix_row_vector) 540 self.lmatrix_col = tuple(lmatrix_col_vector) 541 542 self.lzeromatrix = lomatrix 543 self.lzeromatrix_row = tuple(lomatrix_row_vector) 544 self.lzeromatrix_col = tuple(lomatrix_col_vector) 545 546 self.conservation_matrix = consmatrix 547 self.conservation_matrix_row = tuple(cons_row_vector) 548 self.conservation_matrix_col = tuple(cons_col_vector) 549 550 self.nrmatrix = nrmatrix 551 self.nrmatrix_row = tuple(Nred_vector) 552 self.nrmatrix_col = tuple(scipy.copy(self.nmatrix_col)) 553 print(' done.') 554 555 556 def PivotSort(self,a,row_vector,column_vector): 557 """ 558 PivotSort(a,row_vector,column_vector) 559 560 This is a sorting routine that accepts a matrix and row/colum vectors 561 and then sorts them so that: there are no zero rows (by swapping with first 562 non-zero row) The abs(largest) pivots are moved onto the diagonal to maintain 563 numerical stability. Row and column swaps are recorded in the tracking vectors. 564 565 Arguments: 566 567 a: the input array 568 row_vector: row tracking vector 569 column_vector: column tracking vector 570 571 """ 572 t = self.commonType(a) 573 row,col = a.shape 574 575 for z in range(0,min(row,col)): 576 if abs(a[z,z]) < self.stoichiometric_analysis_lu_precision: 577 maxV = self.stoichiometric_analysis_lu_precision 578 maxP = (None,None) 579 zeroP = (None,None) 580 581 mList = [] 582 for x in range(z,min(row,col)): 583 mList.append((x,max(abs(a[x,x:])))) 584 mVal = 0.0 585 mRow = None 586 #print mList 587 for el in mList: 588 if el[1] > mVal: 589 mVal = el[1] 590 mRow = el[0] 591 if self.__stoichdiagmode__: 592 print('\tmVal', mVal, mRow) 593 594 if mRow != None: 595 for y in range(z,col): 596 if abs(a[mRow,y]) > abs(maxV) and abs(a[mRow,y]) > self.stoichiometric_analysis_lu_precision: 597 maxV = a[mRow,y] 598 maxP = (mRow,y) 599 600 if zeroP != (None,None) and zeroP != (z,z) and mRow != None: 601 if self.__stoichdiagmode__: 602 print(' Swapping: ', (z,z), (zeroP[0],zeroP[1]), 1.0) 603 a = self.SwapRowd(a,z,zeroP[0]) 604 a = self.SwapCold(a,z,zeroP[1]) 605 row_vector = self.SwapElem(row_vector,z,zeroP[0]) 606 column_vector = self.SwapElem(column_vector,z,zeroP[1]) 607 elif maxP != (None,None) and maxP != (min(row,col),min(row,col)) and mRow != None: 608 if self.__stoichdiagmode__: 609 print(' Swapping: ', (z,z), (maxP[0],maxP[1]), a[z,z], maxV) 610 a = self.SwapRowd(a,z,maxP[0]) 611 a = self.SwapCold(a,z,maxP[1]) 612 row_vector = self.SwapElem(row_vector,z,maxP[0]) 613 column_vector = self.SwapElem(column_vector,z,maxP[1]) 614 615 return(a,row_vector,column_vector) 616 617 def PivotSort_initial(self,a,row_vector,column_vector): 618 """ 619 PivotSort_initial(a,row_vector,column_vector) 620 621 This is a sorting routine that accepts a matrix and row/colum vectors 622 and then sorts them so that: the abs(largest) pivots are moved onto the diagonal to maintain 623 numerical stability i.e. the matrix diagonal is in descending max(abs(value)). 624 Row and column swaps are recorded in the tracking vectors. 625 626 Arguments: 627 628 a: the input array 629 row_vector: row tracking vector 630 column_vector: column tracking vector 631 632 """ 633 634 # SAME AS THE ABOVE JUST DOES ALL VALUES NOT ONLY NON_ZERO ONES 635 # brett 200500802 636 637 t = self.commonType(a) 638 row,col = a.shape 639 640 for z in range(0,min(row,col)): 641 if z == 0 or abs(a[z,z]) < abs(a[z-1,z-1]): 642 maxV = self.stoichiometric_analysis_lu_precision 643 maxP = (None,None) 644 zeroP = (None,None) 645 mList = [] 646 for x in range(z,min(row,col)): 647 mList.append((x,max(abs(a[x,x:])))) 648 mVal = 0.0 649 mRow = None 650 for el in mList: 651 if el[1] > mVal: 652 mVal = el[1] 653 mRow = el[0] 654 if self.__stoichdiagmode__: 655 print('\tmVal', mVal, mRow) 656 657 if mRow != None: 658 for y in range(z,col): 659 if abs(a[x,y]) > abs(maxV) and abs(a[x,y]) > self.stoichiometric_analysis_lu_precision: 660 maxV = a[x,y] 661 maxP = (x,y) 662 663 if maxP != (None,None) and maxP != (min(row,col),min(row,col)) and maxP != (z,z): 664 if self.__stoichdiagmode__: 665 print(' Swapping: ', (z,z), (maxP[0],maxP[1]), a[z,z], maxV) 666 a = self.SwapRowd(a,z,maxP[0]) 667 a = self.SwapCold(a,z,maxP[1]) 668 row_vector = self.SwapElem(row_vector,z,maxP[0]) 669 column_vector = self.SwapElem(column_vector,z,maxP[1]) 670 671 return(a,row_vector,column_vector) 672 673 674 def PLUfactorize(self,a_in): 675 """ 676 PLUfactorize(a_in) 677 678 Performs an LU factorization using LAPACK D/ZGetrf. Now optimized for FLAPACK interface. 679 Returns LU - combined factorization, IP - rowswap information and info - Getrf error control. 680 681 Arguments: 682 683 a_in: the matrix to be factorized 684 685 """ 686 print('.', end=' ') 687 688 self.assertRank2(a_in) 689 t = self.commonType(a_in) 690 if a_in.dtype.char == 'D': 691 #print 'Complex matrix ' + a_in.typecode() 692 ## a = copy.copy(scipy.transpose(a_in)) 693 # brett 201106 flapack optimize 694 a = a_in.copy() 695 elif a_in.dtype.char == 'd': 696 #print 'Float matrix ' + a_in.typecode() 697 ## a = copy.copy(scipy.transpose(a_in)) 698 # brett 201106 flapack optimize 699 a = a_in.copy() 700 else: 701 #print 'Other matrix casting to double, was: ' + a_in.typecode() 702 ## a = copy.copy(scipy.transpose(a_in).astype('d')) 703 # brett 201106 flapack optimize 704 a = a_in.copy().astype('d') 705 Using_FLAPACK = 1 706 # brett 20110620 - clapack support is being removed from scipy going exclusively flapack now 707 # brett 20041226 - protecting ourselves against flapack 708 ## try: 709 ## scipy.linalg.clapack.empty_module() 710 ## Using_FLAPACK = 1 711 ## print "Using FLAPACK" 712 ## except: 713 ## print "Using CLAPACK" 714 715 if Using_FLAPACK == 1: 716 try: 717 if self.array_kind[t] == 1: 718 getrf = myflapack.zgetrf #scipy flapack 20041226 719 else: 720 getrf = myflapack.dgetrf #scipy flapack 20041226 721 except Exception as e: 722 print("FLAPACK error", e) 723 ## else: 724 ## try: 725 ## if self.array_kind[t] == 1: 726 ## getrf = scipy.linalg.clapack.zgetrf #scipy clapack (ATLAS) 20030605 727 ## else: 728 ## getrf = scipy.linalg.clapack.dgetrf #scipy clapack (ATLAS) 20030605 729 ## except Exception, e: 730 ## print "CLAPACK error", e 731 732 if Using_FLAPACK == 1: 733 ## results = getrf(scipy.transpose(a)) # brett 20041226 734 results = getrf(a) # brett 201106 735 ## else: 736 ## # This is a $%^&*( ... suddenly with latest cvs f2py getrf only accepts arrays 737 ## # not vectors so this song and dance is necessary to fix this 'behaviour' 738 ## # as idiotic as it seems, we pad the vector with a zero row (no difference numerically) 739 ## # run it through getrf and then strip it afterwards !@#$%^&*() - brett 20050707 740 ## if a.shape[0] == 1: 741 ## tarr = scipy.zeros((a.shape[0]+1,a.shape[1]),'d') 742 ## tarr[0] = a[0] 743 744 ## results = getrf(tarr) #scipy cblas (ATLAS) 20030506 -- this is normal 745 746 ## results = list(results) 747 ## results[0] = scipy.array([results[0][0]]) 748 ## results[1] = scipy.array([results[1][0]],'i') 749 ## results[2] = results[2] 750 ## results = tuple(results) 751 ## else: 752 ## results = getrf(a) #scipy cblas (ATLAS) 20030506 -- this is normal 753 754 results = list(results) 755 756 if results[2] < 0: 757 print('Argument ', results['info'], ' had an illegal value') 758 raise LinAlgError 759 elif results[2] > 0: 760 pass 761 762 if Using_FLAPACK == 1: 763 result = results[0] # brett 20041226 764 ## else: 765 ## result = scipy.transpose(results[0]) # -- this is normal 766 767 badlist = [] 768 for x in range(result.shape[0]): 769 for y in range(result.shape[1]): 770 if abs(result[x,y]) != 0.0 and abs(result[x,y]) < self.stoichiometric_analysis_lu_precision: 771 result[x,y] = 0.0 # gets rid of the -0.0 irritation 772 if len(badlist) == 0 and (x,y) == (x,x): # catch 1st float on a pivot 773 badlist.append(x) 774 if len(badlist) != 0: 775 results[2] = badlist[0]-1 # 20040423 if float was on a pivot - refactorize 776 777 return(result,results[1],results[2]) #scipy cblas (ATLAS?) 20030506 778 779 780 ## def PLUfactorizeOLD(self,a_in): 781 ## """ 782 ## PLUfactorize(a_in) OLD PRE SCIPY 0.10 uses clapack 783 784 ## Performs an LU factorization using LAPACK D/ZGetrf. 785 ## Returns LU - combined factorization, IP - rowswap information and info - Getrf error control. 786 787 ## Arguments: 788 789 ## a_in: the matrix to be factorized 790 791 ## """ 792 ## print '.', 793 794 ## self.assertRank2(a_in) 795 ## t = self.commonType(a_in) 796 ## if a_in.dtype.char == 'D': 797 ## #print 'Complex matrix ' + a_in.typecode() 798 ## a = copy.copy(scipy.transpose(a_in)) 799 ## elif a_in.dtype.char == 'd': 800 ## #print 'Float matrix ' + a_in.typecode() 801 ## a = copy.copy(scipy.transpose(a_in)) 802 ## else: 803 ## #print 'Other matrix casting to double, was: ' + a_in.typecode() 804 ## a = copy.copy(scipy.transpose(a_in).astype('d')) 805 ## Using_FLAPACK = 0 806 ## # brett 20041226 - protecting ourselves against flapack 807 ## try: 808 ## scipy.linalg.clapack.empty_module() 809 ## Using_FLAPACK = 1 810 ## print "Using FLAPACK" 811 ## except: 812 ## print "Using CLAPACK" 813 814 ## if Using_FLAPACK == 1: 815 ## try: 816 ## if self.array_kind[t] == 1: 817 ## getrf = myflapack.zgetrf #scipy flapack 20041226 818 ## else: 819 ## getrf = myflapack.dgetrf #scipy flapack 20041226 820 ## except Exception, e: 821 ## print "FLAPACK error", e 822 ## else: 823 ## try: 824 ## if self.array_kind[t] == 1: 825 ## getrf = scipy.linalg.clapack.zgetrf #scipy clapack (ATLAS) 20030605 826 ## else: 827 ## getrf = scipy.linalg.clapack.dgetrf #scipy clapack (ATLAS) 20030605 828 ## except Exception, e: 829 ## print "CLAPACK error", e 830 831 ## if Using_FLAPACK == 1: 832 ## results = getrf(scipy.transpose(a)) # brett 20041226 833 ## else: 834 ## # This is a $%^&*( ... suddenly with latest cvs f2py getrf only accepts arrays 835 ## # not vectors so this song and dance is necessary to fix this 'behaviour' 836 ## # as idiotic as it seems, we pad the vector with a zero row (no difference numerically) 837 ## # run it through getrf and then strip it afterwards !@#$%^&*() - brett 20050707 838 ## if a.shape[0] == 1: 839 ## tarr = scipy.zeros((a.shape[0]+1,a.shape[1]),'d') 840 ## tarr[0] = a[0] 841 842 ## results = getrf(tarr) #scipy cblas (ATLAS) 20030506 -- this is normal 843 844 ## results = list(results) 845 ## results[0] = scipy.array([results[0][0]]) 846 ## results[1] = scipy.array([results[1][0]],'i') 847 ## results[2] = results[2] 848 ## results = tuple(results) 849 ## else: 850 ## results = getrf(a) #scipy cblas (ATLAS) 20030506 -- this is normal 851 852 ## results = list(results) 853 854 ## if results[2] < 0: 855 ## print('Argument ', results['info'], ' had an illegal value') 856 ## raise LinAlgError 857 ## elif results[2] > 0: 858 ## pass 859 860 ## if Using_FLAPACK == 1: 861 ## result = results[0] # brett 20041226 862 ## else: 863 ## result = scipy.transpose(results[0]) # -- this is normal 864 865 866 ## badlist = [] 867 ## for x in range(result.shape[0]): 868 ## for y in range(result.shape[1]): 869 ## if abs(result[x,y]) != 0.0 and abs(result[x,y]) < self.stoichiometric_analysis_lu_precision: 870 ## result[x,y] = 0.0 # gets rid of the -0.0 irritation 871 ## if len(badlist) == 0 and (x,y) == (x,x): # catch 1st float on a pivot 872 ## badlist.append(x) 873 ## if len(badlist) != 0: 874 ## results[2] = badlist[0]-1 # 20040423 if float was on a pivot - refactorize 875 876 ## return(result,results[1],results[2]) #scipy cblas (ATLAS?) 20030506 877 878 def SplitLU(self,plu,row,col,t=None): 879 """ 880 SplitLU(plu,row,col,t) 881 882 PLU takes the combined LU factorization computed by PLUfactorize and extracts the upper matrix. 883 Returns U. 884 885 Arguments: 886 887 plu: LU factorization 888 row: row tracking vector 889 col: column tracking vector 890 t [default=None)]: typecode argument (currently not used) 891 892 """ 893 print('.', end=' ') 894 895 if self.__stoichdiagmode__: 896 print('\n',row,col) 897 for cl in range(0,min(row,col)): 898 plu[cl+1:,cl] = 0.0 899 return plu 900 901 def GetUpperMatrix(self,a): 902 """ 903 GetUpperMatrix(a) 904 905 Core analysis algorithm; an input is preconditioned using PivotSort_initial and then cycles of PLUfactorize and 906 PivotSort are run until the factorization is completed. During this process the matrix is reordered by 907 column swaps which emulates a full pivoting LU factorization. Returns the pivot matrix P, upper factorization U 908 as well as the row/col tracking vectors. 909 910 Arguments: 911 912 a: a stoichiometric matrix 913 914 """ 915 print('.', end=' ') 916 917 t = self.commonType(a) 918 row,col = a.shape 919 920 # this is a test brett 20050802 921 row_vector = scipy.array((list(range(row)))) 922 column_vector = scipy.array((list(range(col)))) 923 a,row_vector,column_vector = self.PivotSort_initial(a,row_vector,column_vector) 924 925 #18/09/2000 PLUfactorize included in self.GetUpperMatrix 926 a,ip,info = self.PLUfactorize(a) 927 928 # get U from getrf's PLU 929 if self.__stoichdiagmode__: 930 print(info) 931 upper_out = self.SplitLU(a,row,col,t) 932 933 p_out = scipy.identity(row) 934 for x in range(0,min(row,col)): 935 if x + 1 != ip[x]: 936 p_out = self.SwapRowd(p_out,x,(ip[x]-1)) 937 row_vector = self.SwapElem(row_vector,x,(ip[x]-1)) 938 939 '''31/08/2000 Added to try to make sure that pivots exist only on the upper left side of the matrix 940 max(abs(upper_out[x,:])) used instead of abs(max(upper_out[x,:])) otherwise for rows where 941 all elements < 0 zero is maximum, this way the max of positive values is used''' 942 943 upper_out,row_vector,column_vector = self.PivotSort(upper_out,row_vector,column_vector) 944 945 946 '''20/09/2000 This bit sorts out the echelon matrix by running (if needed) cycles of 947 self.PLUfactorize, self.GetUpperMatrix, and pivsort until only a staircase matrix remains. It uses both the error 948 generated by self.PLUfactorize(info) and go_flag to control itself''' 949 950 if info > 0: 951 #Here we check to see if the error is not in the last or reduced-last column if it is - exit 952 # go_flag = 'yes' # 2001/04/26 changed for Python21 and future compatibility 953 go_flag = 1 954 for x in range (0,min(row,col)): 955 if abs(upper_out[x,x]) > self.stoichiometric_analysis_lu_precision: 956 pos_holder = x 957 if pos_holder + 1 < info and pos_holder + 1 < row: 958 if max(abs(upper_out[pos_holder + 1,:])) < self.stoichiometric_analysis_lu_precision: 959 # go_flag = 'no' # 2001/04/26 changed for Python21 and future compatibility 960 go_flag = 0 961 elif pos_holder + 1 == info and pos_holder + 1 == min(row,col): 962 # go_flag = 'no' # 2001/04/26 changed for Python21 and future compatibility 963 go_flag = 0 964 # while info > 0 and go_flag == 'yes': # 2001/04/26 changed for Python21 and future compatibility 965 while info > 0 and go_flag == 1: 966 #sort 967 #upper_out,row_vector,column_vector = pivot_sort2k5(upper_out,row_vector,column_vector) 968 upper_out,row_vector,column_vector = self.PivotSort(upper_out,row_vector,column_vector) 969 #PLUfactorize 970 if self.__stoichdiagmode__: 971 print(' Echelon: ' + time.strftime(self.__TimeFormat), '(', info,')') 972 upper_out,ip,info = self.PLUfactorize(upper_out) 973 974 # get U from getrf's PLU 975 if self.__stoichdiagmode__: 976 print(info) 977 upper_out = self.SplitLU(upper_out,row,col,t) 978 979 #realign row vector and permutation matrix if necessary 980 for x in range(0,min(row,col)): 981 if x + 1 != ip[x]: 982 p_out = self.SwapRowd(p_out,x,(ip[x]-1)) 983 row_vector = self.SwapElem(row_vector,x,(ip[x]-1)) 984 #test if we can exit -- same criteria as before 985 for x in range (0,min(row,col)): 986 if abs(upper_out[x,x]) > self.stoichiometric_analysis_lu_precision: 987 pos_holder = x 988 if pos_holder + 1 < info and pos_holder + 1 < row: 989 if max(abs(upper_out[pos_holder + 1,:])) < self.stoichiometric_analysis_lu_precision: 990 # go_flag = 'no' # 2001/04/26 changed for Python21 and future compatibility 991 go_flag = 0 992 elif pos_holder + 1 == info and pos_holder + 1 == min(row,col): 993 # go_flag = 'no' # 2001/04/26 changed for Python21 and future compatibility 994 go_flag = 0 995 996 '''This bit will get rid of any zero rows so that we will hopefully only have to work with a 997 reduced matrix after this, for completeness the row_vector will also be sliced''' 998 999 for x in range (0,min(row,col)): 1000 if abs(upper_out[x,x]) > self.stoichiometric_analysis_lu_precision: 1001 pos_holder = x 1002 1003 upper_out_r = scipy.zeros((pos_holder + 1,col)).astype(t) 1004 upper_out_r = upper_out[:pos_holder + 1,:] 1005 row_vector_r = row_vector[:pos_holder + 1] 1006 1007 return(p_out,upper_out_r,row_vector_r,column_vector) 1008 1009 1010 def GetUpperMatrixUsingQR(self,a): 1011 """ 1012 GetUpperMatrix(a) 1013 1014 Core analysis algorithm; an input is preconditioned using PivotSort_initial and then cycles of PLUfactorize and 1015 PivotSort are run until the factorization is completed. During this process the matrix is reordered by 1016 column swaps which emulates a full pivoting LU factorization. Returns the pivot matrix P, upper factorization U 1017 as well as the row/col tracking vectors. 1018 1019 Arguments: 1020 1021 a: a stoichiometric matrix 1022 1023 """ 1024 print('.', end=' ') 1025 1026 t = self.commonType(a) 1027 row,col = a.shape 1028 1029 # this is a test brett 20050802 1030 row_vector = scipy.array((list(range(row)))) 1031 column_vector = scipy.array((list(range(col)))) 1032 ## a,row_vector,column_vector = self.PivotSort(a,row_vector,column_vector) 1033 a,row_vector,column_vector = self.PivotSort_initial(a,row_vector,column_vector) 1034 1035 1036 Q,upper_out = scipy.linalg.qr(a) 1037 self.MatrixFloatFix(upper_out,val=self.stoichiometric_analysis_lu_precision*10.0) 1038 1039 1040 '''This bit will get rid of any zero rows so that we will hopefully only have to work with a 1041 reduced matrix after this, for completeness the row_vector will also be sliced''' 1042 1043 for x in range (0,min(row,col)): 1044 if abs(upper_out[x,x]) > self.stoichiometric_analysis_lu_precision: 1045 pos_holder = x 1046 1047 1048 upper_out_r = scipy.zeros((pos_holder + 1,col)).astype(t) 1049 p_out = scipy.diag(upper_out.shape[1]*[1.0]) 1050 upper_out_r = upper_out[:pos_holder + 1,:] 1051 row_vector_r = row_vector[:pos_holder + 1] 1052 1053 return(p_out,upper_out_r,row_vector_r,column_vector) 1054 1055 def ScalePivots(self,a_one): 1056 """ 1057 ScalePivots(a_one) 1058 1059 Given an upper triangular matrix U, this method scales the diagonal (pivot values) to one. 1060 1061 Arguments: 1062 1063 a_one: an upper triangular matrix U 1064 1065 """ 1066 print('.', end=' ') 1067 1068 t = self.commonType(a_one) 1069 row,col = a_one.shape 1070 1071 '''13/09/2000 We now assume that the matrix has the correct shape ie. the pivots are in a 1072 perfect staircase''' 1073 1074 for x in range (0,row): 1075 if abs(a_one[x,x]) == 0.0: # ths is to prevent NAN's when FixFloat zeros a supersmall pivot 1076 pass 1077 elif abs(a_one[x,x]) != 1.0: 1078 a_one[x,:] = a_one[x,:]/abs(a_one[x,x]) 1079 if a_one[x,x] < 0.0: 1080 a_one[x,:] = -(a_one[x,:]) 1081 return(a_one) 1082 1083 1084 def BackSubstitution(self,res_a,row_vector,column_vector): 1085 """ 1086 BackSubstitution(res_a,row_vector,column_vector) 1087 1088 Jordan reduction of a scaled upper triangular matrix. The returned array is now in the form [I R] and can 1089 be used for nullspace determination. Modified row and column tracking vetors are also returned. 1090 1091 Arguments: 1092 1093 res_a: unitary pivot upper triangular matrix 1094 row_vector: row tracking vector 1095 column_vector: column tracking vector 1096 1097 """ 1098 print('.', end=' ') 1099 1100 t = self.commonType(res_a) 1101 row,col = res_a.shape 1102 1103 '''13/09/2000 removed the copy thing, and eliminated the divnumb variable. Because we 1104 are now working with a row reduced matrix upward elimination only works on the pivot cols''' 1105 1106 # old back substitution circa 2000! brett - 20050801 1107 ## bigF = 0 1108 ## if max(res_a.shape) > 500: 1109 ## bigF = 1 1110 ## for y in range (min(row,col)-1,-1,-1): 1111 ## if bigF and self.__stoichdiagmode__: 1112 ## print '.', 1113 ## for x in range (min(row,col)-2,-1,-1): 1114 ## z = x + 1 1115 ## while z < min(row,col): 1116 ## if abs(res_a[x,y]) > self.stoichiometric_analysis_gj_precision and abs(res_a[z,y]) > self.stoichiometric_analysis_gj_precision: 1117 ## res_a[x,:] = res_a[x,:] - res_a[z,:]*(res_a[x,y]/res_a[z,y]) 1118 ## z = z + 1 1119 ## continue 1120 ## else: 1121 ## z = z + 1 1122 1123 # new back substitution 1124 # right looking algorithm that uses array slicing speeds up as it moves down the pivots 1125 # this algorithm is at least 3X as fast as the old one for a large system - brett 20050805 1126 bigF = 0 1127 if max(res_a.shape) > 500: 1128 bigF = 1 1129 for x in range(min(row,col)): 1130 if bigF and self.__stoichdiagmode__: 1131 print('.', end=' ') 1132 for y in range(x+1,min(row,col)): 1133 if abs(res_a[x,y]) > self.stoichiometric_analysis_lu_precision: 1134 res_a[x,y:] = res_a[x,y:]-(res_a[x,y]*res_a[y,y:]) 1135 1136 # $%^&* wtf is this for??? brett - 20050801 (besides cleaning fp wierdness ... nothing ;-) 1137 self.MatrixFloatFix(res_a,val=self.stoichiometric_analysis_gj_precision) 1138 res_a,row_vector,column_vector = self.PivotSort(res_a,row_vector,column_vector) 1139 1140 return (res_a,row_vector,column_vector) 1141 1142 1143 def K_split_R(self,R_a,row_vector,column_vector): 1144 """ 1145 K_split_R(R_a,row_vector,column_vector) 1146 1147 Using the R factorized form of the stoichiometric matrix we now form the K and Ko matrices. Returns 1148 the r_ipart,Komatrix,Krow,Kcolumn,Kmatrix,Korow,info 1149 1150 Arguments: 1151 1152 R_a: the Gauss-Jordan reduced stoichiometric matrix 1153 row_vector: row tracking vector 1154 column_vector: column tracking vector 1155 1156 """ 1157 print('.', end=' ') 1158 1159 t = self.commonType(R_a) 1160 row,col = R_a.shape 1161 1162 '''14/09/2000 Seeing as we now should have a perfect staircase in a reduced matrix we do 1163 not have to search for the last pivot it will exist at min(row,col)-1 so the pos_holder 1164 finding code has been removed (actually moved into self.GetUpperMatrix)''' 1165 1166 pos_holder = min(row,col) - 1 1167 1168 '''This bit extracts the identity part from R (future note this could be replaced by an I matrix formed by min(row,col))''' 1169 1170 r_ipart = scipy.zeros((pos_holder + 1,pos_holder + 1)).astype(t) 1171 r_ipart = R_a[:pos_holder+1,:pos_holder+1] 1172 1173 row_i,col_i = r_ipart.shape 1174 1175 '''If there are free variables, then this bit will extract them and form the row/col vectors''' 1176 1177 empty_rf = 0 1178 K_switch = 0 #added 20020416 class attribute for conservation detection 0 = none, 1 = exists 1179 1180 if col - col_i > self.stoichiometric_analysis_fp_zero: 1181 r_fpart = scipy.zeros((pos_holder + 1,col - (pos_holder + 1))).astype(t) 1182 r_fpart = R_a[:pos_holder+1,pos_holder+1:] 1183 1184 row_vector_dependent = column_vector[:pos_holder + 1] 1185 row_vector_independent = column_vector[pos_holder + 1:] 1186 1187 #row_vector = scipy.concatenate((row_vector_independent,row_vector_dependent),1) 1188 row_vector = scipy.hstack((row_vector_independent,row_vector_dependent)) # numpy 0.10 1189 column_vector = column_vector[pos_holder + 1:] 1190 K_switch = 1 1191 else: 1192 print('no flux conservation') 1193 r_fpart = 'no flux conservation' 1194 empty_rf = 1 1195 K_switch = 0 1196 1197 # if r_fpart == 'F is empty, R = I': # 2001/04/26 changed for Python21 compatibility 1198 if empty_rf == 1: 1199 return(r_ipart,r_ipart,column_vector,column_vector,r_ipart,column_vector,K_switch) 1200 else: 1201 row,col = r_fpart.shape 1202 id = scipy.identity(col).astype(t) 1203 nullspace = scipy.concatenate((id,-r_fpart),0).astype(t) 1204 1205 #brett 05/11/2002 changed r_fpart to -r_fpart 1206 return(r_ipart,-r_fpart,row_vector,column_vector,nullspace,row_vector_dependent,K_switch) 1207 1208 1209 def L_split_R(self,Nfull,R_a,row_vector,column_vector): 1210 """ 1211 L_split_R(Nfull,R_a,row_vector,column_vector) 1212 1213 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 1214 1215 Arguments: 1216 1217 Nfull: the original stoichiometric matrix N 1218 R_a: gauss-jordan factorized form of N^T 1219 row_vector: row tracking vector 1220 column_vector: column tracking vector 1221 1222 """ 1223 print('.', end=' ') 1224 1225 #print '\nR_a' 1226 #print `R_a` 1227 t = self.commonType(R_a) 1228 row,col = R_a.shape 1229 1230 '''14/09/2000 Seeing as we now should have a perfect staircase in a reduced matrix we do 1231 not have to search for the last pivot it will exist at min(row,col)-1 so the pos_holder 1232 finding code has been removed (actually moved into self.GetUpperMatrix)''' 1233 1234 pos_holder = min(row,col) - 1 1235 1236 '''Here we extract the identity matrix from R (future note this could be replaced by an I matrix formed by min(row,col))''' 1237 1238 r_ipart = scipy.zeros((pos_holder + 1,pos_holder + 1)).astype(t) 1239 r_ipart = R_a[:pos_holder+1,:pos_holder+1] 1240 1241 row_i,col_i = r_ipart.shape 1242 1243 # exit1 = 'no' # 2001/04/26 changed for Python21 and future compatibility 1244 exit1 = 0 1245 L_switch = False #added 20020416 class attribute for conservation detection 0 = none, 1 = exists 1246 1247 '''If there are free variable then they are extracted and packaged, row/col vectors are formed''' 1248 1249 if col - col_i > self.stoichiometric_analysis_fp_zero: 1250 r_fpart = scipy.zeros((pos_holder + 1,col - (pos_holder + 1))).astype(t) 1251 r_fpart = R_a[:pos_holder+1,pos_holder+1:] 1252 1253 col_vector_dependent = column_vector[pos_holder + 1:] 1254 col_vector_independent = column_vector[:pos_holder + 1] 1255 #cons_col_vector = scipy.concatenate((col_vector_independent,col_vector_dependent),1) 1256 cons_col_vector = scipy.hstack((col_vector_independent,col_vector_dependent)) # numpy 0.10 1257 cons_row_vector = col_vector_dependent 1258 #lmatrix_row_vector = scipy.concatenate((col_vector_independent,col_vector_dependent),1) 1259 lmatrix_row_vector = scipy.hstack((col_vector_independent,col_vector_dependent)) # numpy 0.10 1260 lmatrix_col_vector = col_vector_independent 1261 lomatrix_row_vector = col_vector_dependent 1262 lomatrix_col_vector = col_vector_independent 1263 Nred_vector = col_vector_independent 1264 L_switch = True 1265 else: 1266 lomatrix_row_vector = column_vector[:pos_holder + 1] 1267 lomatrix_col_vector = column_vector[:pos_holder + 1] 1268 Nred_vector = column_vector[:pos_holder + 1] 1269 #print('F is empty, R = I, no conservation matrix, L is Lo') 1270 r_fpart = r_ipart 1271 # exit1 = 'yes' # 2001/04/26 changed for Python21 and future compatibility 1272 exit1 = 1 1273 L_switch = False 1274 1275 # if exit1 == 'yes': # 2001/04/26 changed for Python21 and future compatibility 1276 if exit1 == 1: 1277 r_fpart = scipy.transpose(r_fpart) 1278 lmatrix = r_fpart 1279 #02/10/2000 removed so that the thing returns the Lo matrix as L 1280 #r_fpart = 'no conservation Lo = L = I' 1281 # my factorization routines now swap things around for numeric stability, so that Nr might be a row/column swapped 1282 # this simply synchronizes Nr with its labels - brett 20050805 1283 Nfull = scipy.transpose(Nfull) 1284 row,col = Nfull.shape 1285 Nred = scipy.zeros((row_i,col)).astype(t) 1286 for x in range (0,row_i): 1287 Nred[x,:] = Nfull[Nred_vector[x],:] 1288 # return the right stuff - brett 20050805 1289 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) 1290 #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) 1291 else: 1292 r_fpart = scipy.transpose(r_fpart) 1293 row,col = r_fpart.shape 1294 1295 id = scipy.identity(row).astype(t) 1296 #consmatrix = scipy.concatenate((-r_fpart,id),1).astype(t) 1297 consmatrix = scipy.hstack((-r_fpart,id)).astype(t) # numpy 0.10 1298 1299 id = scipy.identity(col).astype(t) 1300 lmatrix = scipy.concatenate((id,r_fpart),0).astype(t) 1301 1302 '''This bit creates Nr. The transpose is only necessary if Nfull is already transposed in the input function''' 1303 1304 Nfull = scipy.transpose(Nfull) 1305 row,col = Nfull.shape 1306 Nred = scipy.zeros((row_i,col)).astype(t) 1307 for x in range (0,row_i): 1308 Nred[x,:] = Nfull[Nred_vector[x],:] 1309 1310 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) 1311 1312 1313 def SVD_Rank_Check(self,matrix=None,factor=1.0e4,resultback=0): 1314 """ 1315 SVD_Rank_Check(matrix=None,factor=1.0e4,resultback=0) 1316 1317 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 1318 1319 Arguments: 1320 1321 matrix [default=None]: the stoichiometric matrix default is self.Nmatrix 1322 factor [default=1.0e4]: factor used to calculate the 'zero pivot' mask = mach_eps*factor 1323 resultback [default=0]: return the SVD results, U, S, vh 1324 1325 """ 1326 if matrix == None: 1327 matrix = self.nmatrix 1328 1329 nrow,ncol = matrix.shape 1330 TrMat = 0 1331 if ncol > nrow: 1332 # 'INFO: SVD is more accurate for tall matrices using (N)T\n' 1333 matrix = scipy.transpose(matrix) 1334 TrMat = 1 1335 1336 u,s,vh = scipy.linalg.svd(matrix) 1337 maskF = scipy.machar.machar_double.eps*factor 1338 1339 if TrMat == 0: 1340 print('SVD zero mask:', maskF) 1341 print('LU factorization effective zero:', self.stoichiometric_analysis_lu_precision) 1342 1343 rank = len([el for el in (abs(s) > maskF) if el > 0.0]) 1344## null_mask = (abs(s) > maskF) 1345## vhnull = scipy.compress(null_mask,vh,axis=0) 1346## unull = scipy.compress(null_mask,u,axis=0) 1347## assert vhnull.shape[0] == unull.shape[0], 'This should be true or I\'m very confused' 1348## rank = vhnull.shape[0] 1349 1350 1351 print('\nNmatrix has ' + repr(nrow) + ' rows and ' + repr(ncol) + ' columns') 1352 1353 print('\nSVD \"considers\" the rank to be: ' + repr(rank)) 1354 print('LU (Kmatrix) considers the rank to be: ' + repr(self.kzeromatrix.shape[0])) 1355 print('LU (Lmatrix) considers the rank to be: ' + repr(self.lzeromatrix.shape[1])) 1356 1357 print('\nComparing lzeromatrix dimensions') 1358 print('SVD lzeromatrix.shape =', (u.shape[0]-rank,rank)) 1359 print('LU lzeromatrix.shape =', self.lzeromatrix.shape) 1360 1361 print('\nComparing lmatrix dimensions') 1362 print('SVD lmatrix.shape =', (u.shape[0],rank)) 1363 print('LU lmatrix.shape =', self.lmatrix.shape) 1364 1365 print('\nComparing kzeromatrix dimensions') 1366 print('SVD kzeromatrix.shape =', (rank,vh.shape[0]-rank)) 1367 print('LU kzeromatrix.shape =', self.kzeromatrix.shape) 1368 1369 print('\nComparing kmatrix dimensions') 1370 print('SVD kmatrix.shape =', (vh.shape[0],vh.shape[0]-rank)) 1371 print('LU kmatrix.shape =', self.kmatrix.shape) 1372 else: 1373 print('SVD zero mask:', maskF) 1374 print('LU factorization effective zero:', self.stoichiometric_analysis_lu_precision) 1375 1376 rank = len([el for el in (abs(s) > maskF) if el > 0.0]) 1377## null_mask = (abs(s) > maskF) 1378## vhnull = scipy.compress(null_mask,vh,axis=0) 1379## unull = scipy.compress(null_mask,u,axis=0) 1380## assert vhnull.shape[0] == unull.shape[0], 'This should be true or I\'m very confused' 1381## rank = vhnull.shape[0] 1382 1383 1384 print('\nNmatrix has ' + repr(nrow) + ' rows and ' + repr(ncol) + ' columns') 1385 1386 print('\nSVD considers the rank to be: ' + repr(rank)) 1387 print('LU (Kmatrix) considers the rank to be: ' + repr(self.kzeromatrix.shape[0])) 1388 print('LU (Lmatrix) considers the rank to be: ' + repr(self.lzeromatrix.shape[1])) 1389 1390 print('\nComparing lzeromatrix dimensions') 1391 print('SVD lzeromatrix.shape =', (vh.shape[0]-rank,rank)) 1392 print('LU lzeromatrix.shape =', self.lzeromatrix.shape) 1393 1394 print('\nComparing lmatrix dimensions') 1395 print('SVD lmatrix.shape =', (vh.shape[0],rank)) 1396 print('LU lmatrix.shape =', self.lmatrix.shape) 1397 1398 print('\nComparing kzeromatrix dimensions') 1399 print('SVD kzeromatrix.shape =', (rank,u.shape[0]-rank)) 1400 print('LU kzeromatrix.shape =', self.kzeromatrix.shape) 1401 1402 print('\nComparing kmatrix dimensions') 1403 print('SVD kmatrix.shape =', (u.shape[0],u.shape[0]-rank)) 1404 print('LU kmatrix.shape =', self.kmatrix.shape) 1405 1406 print('\nMatrix value check\n******************') 1407 1408 print('\nKzeromatrix:') 1409 sm,bg = self.MatrixValueCompare(self.kzeromatrix) 1410 print('Smallest value: ', abs(sm)) 1411 print('Largest value: ', abs(bg)) 1412 print('Ratio abs(bg/sm):', abs(bg/sm)) 1413 1414 print('\nLzeromatrix:') 1415 sm,bg = self.MatrixValueCompare(self.lzeromatrix) 1416 print('Smallest value: ', abs(sm)) 1417 print('Largest value: ', abs(bg)) 1418 print('Ratio abs(bg/sm):', abs(bg/sm)) 1419 1420 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.') 1421 print(' ') 1422 1423 if resultback: 1424 return u,s,vh 1425 1426 1427if __psyco_active__: 1428 psyco.bind(MathArrayFunc) 1429 psyco.bind(Stoich) 1430 1431