1""" 2PySCeS - Python Simulator for Cellular Systems (http://pysces.sourceforge.net) 3 4Copyright (C) 2004-2022 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""" 17 18from __future__ import division, print_function 19from __future__ import absolute_import 20from __future__ import unicode_literals 21 22from .version import __version__ 23 24import os 25import math, operator 26import numpy 27 28from .. import output_dir, __CHGDIR_ON_START__ 29if __CHGDIR_ON_START__: 30 CWD = output_dir 31else: 32 CWD = os.getcwd() 33 34from .InfixParser import MyInfixParser 35 36InfixParser = MyInfixParser() 37InfixParser.buildlexer() 38InfixParser.buildparser( 39 debug=0, 40 debugfile='infix.dbg', 41 tabmodule='infix_tabmodule', 42 outputdir=output_dir, 43) 44InfixParser.setNameStr('self.', '()') 45 46os.chdir(CWD) 47 48class MapList(list): 49 def __init__(self, *args): 50 list.__init__(self, *args) 51 52 def asSet(self): 53 return set(self.__getslice__(0, self.__len__())) 54 55 56class NewCoreBase(object): 57 __DEBUG__ = False 58 name = None 59 annotations = None 60 61 def getName(self): 62 return self.name 63 64 def setName(self, name): 65 self.name = name 66 67 def get(self, attr): 68 """Return an attribute whose name is str(attr)""" 69 return self.__getattribute__(attr) 70 71 def getAnnotation(self): 72 """Returns an annotation dictionary""" 73 if self.annotations == None: 74 self.annotations = {} 75 return self.annotations.copy() 76 77 def setAnnotation(self, key, value): 78 """Set an annotation as a key:value pair""" 79 if self.annotations == None: 80 self.annotations = {} 81 self.annotations.update({key: value}) 82 83 84class NumberBase(NewCoreBase): 85 value = None 86 value_initial = None 87 88 def __call__(self): 89 return self.value 90 91 def getValue(self): 92 return self.value 93 94 def setValue(self, v): 95 self.value = v 96 97 98class Compartment(NewCoreBase): 99 size = None 100 dimensions = None 101 Compartment = None 102 reactions = None 103 species = None 104 area = None 105 106 def __init__(self, name, compartment=None): 107 self.name = name 108 self.Compartment = compartment 109 self.reactions = [] 110 self.species = [] 111 112 def __call__(self): 113 return self.size 114 115 def setSize(self, size, dim): 116 self.size = size 117 assert dim in [0, 1, 2, 3], '\nOkeee! %s dimensions?' % dim 118 self.dimensions = dim 119 120 def setArea(self, area=None): 121 if area == None and self.dimensions == 2: 122 self.area = self.size 123 if self.__DEBUG__: 124 print('Setting reactive area to size for 2D compartment %s' % self.name) 125 elif area == None and self.dimensions == 3: 126 self.area = (113.09733552923255 * self.size ** 2.0) ** (0.33333333333333331) 127 if self.__DEBUG__: 128 print( 129 'Setting reactive area to surface area for 3D compartment %s (assuming a sphere geometry)' 130 % self.name 131 ) 132 self.area = area 133 134 def hasReactions(self): 135 return MapList([r.name for r in self.reactions]) 136 137 def hasSpecies(self): 138 return MapList([s.name for s in self.species]) 139 140 def addReaction(self, reaction): 141 if reaction.name not in self.hasReactions(): 142 self.reactions.append(reaction) 143 self.__setattr__(reaction.name, reaction) 144 if self.__DEBUG__: 145 print('Adding reaction %s' % reaction.name) 146 147 def addSpecies(self, species): 148 if species.name not in self.hasSpecies(): 149 self.species.append(species) 150 self.__setattr__(species.name, species) 151 if self.__DEBUG__: 152 print('Adding species %s' % species.name) 153 else: 154 if self.__DEBUG__: 155 print('Species %s already added' % species.name) 156 157 def getDimensions(self): 158 return self.dimensions 159 160 def getCompartment(self): 161 return self.Compartment 162 163 def hasCompartment(self): 164 if self.Compartment != None: 165 return True 166 else: 167 return False 168 169 def isVolume(self): 170 if self.dimensions == 3: 171 return True 172 else: 173 return False 174 175 def isArea(self): 176 if self.dimensions == 2: 177 return True 178 else: 179 return False 180 181 def isLength(self): 182 if self.dimensions == 1: 183 return True 184 else: 185 return False 186 187 def isPoint(self): 188 if self.dimensions == 0: 189 return True 190 else: 191 return False 192 193 194class BaseUnit(NewCoreBase): 195 '''Base Unit can be of type: time, substance, volume''' 196 197 _types = ('time', 'substance', 'volume', 'area', 'length') 198 value = 1.0 199 type = None 200 201 def __init__(self, name, type): 202 self.name = name 203 assert type in self._types, '\nType must be one of: %s' % str(self._types) 204 self.type = type 205 206 def __call__(self): 207 return self.value 208 209 def getType(self): 210 return self.type 211 212 213class SimpleUnit(NewCoreBase): 214 exponent = 1.0 215 scale = 0.0 216 multiplier = 1.0 217 baseunit = None 218 type = None 219 220 def __init__(self, baseunit, name, exp=1.0, scale=0.0, mult=1.0): 221 self.baseunit = baseunit 222 self.exponent = exp 223 self.scale = scale 224 self.multiplier = mult 225 self.name = name 226 self.type = baseunit.type 227 228 def __call__(self): 229 return (self.multiplier * self.baseunit() * 10 ** self.scale) ** self.exponent 230 231 def getType(self): 232 return self.type 233 234 235class CompoundUnit(NewCoreBase): 236 units = None 237 _HAS_USERNAME = False 238 239 def __init__(self, name=None): 240 self.units = [] 241 if name != None: 242 self.name = name 243 self._HAS_USERNAME = True 244 else: 245 self.name = '' 246 247 def __call__(self): 248 U = 1.0 249 for u in self.units: 250 U *= u() 251 return U 252 253 def addUnit(self, unit): 254 self.units.append(unit) 255 if not self._HAS_USERNAME: 256 self.name = '%s%s' % (self.name, unit.getName()) 257 258 def getUnits(self): 259 return self.units 260 261 def hasUnits(self): 262 return MapList([u.getName() for u in self.units]) 263 264 265class Species(NumberBase): 266 subs = None 267 prods = None 268 mods = None 269 fixed = False 270 Compartment = None 271 __amount__ = False 272 273 def __init__(self, name, value): 274 self.setName(name) 275 self.value = value 276 self.value_initial = value 277 self.subs = [] 278 self.prods = [] 279 self.mods = [] 280 281 def getCompartment(self): 282 return self.Compartment 283 284 def setCompartment(self, c): 285 self.Compartment = c 286 287 def hasCompartment(self): 288 if self.Compartment != None: 289 return True 290 else: 291 return False 292 293 def setSubstrate(self, reaction): 294 self.__setattr__(reaction.name, reaction) 295 self.subs.append(reaction) 296 297 def setProduct(self, reaction): 298 self.__setattr__(reaction.name, reaction) 299 self.prods.append(reaction) 300 301 def setModifier(self, reaction): 302 self.__setattr__(reaction.name, reaction) 303 self.mods.append(reaction) 304 305 def isSubstrateOf(self): 306 return MapList([r.name for r in self.subs]) 307 308 def isProductOf(self): 309 return MapList([r.name for r in self.prods]) 310 311 def isModifierOf(self): 312 return MapList([r.name for r in self.mods]) 313 314 def isReagentOf(self): 315 return MapList(self.isSubstrateOf() + self.isProductOf()) 316 317 def setAmount(self, b): 318 self.__amount__ = bool(b) 319 320 def isAmount(self): 321 return self.__amount__ 322 323 324class SpeciesAssignmentRule(Species): 325 formula = None 326 code_string = None 327 _names = None 328 _functions = None 329 type = 'assignment' 330 _TIME_ = None 331 332 def __init__(self, name, value): 333 Species.__init__(self, name, value) 334 335 def __call__(self): 336 exec(self.xcode) 337 return self.value 338 339 def addFormula(self, formula): 340 formula = formula.replace('self.', '') 341 self.formula = formula 342 InfixParser.setNameStr('self.', '()') 343 InfixParser.parse(formula) 344 self.code_string = 'self.value=%s' % InfixParser.output 345 self._names = InfixParser.names 346 self._functions = InfixParser.functions 347 self.xcode = compile(self.code_string, '<string>', 'exec') 348 349 def addModelAttr(self, obj): 350 self.__setattr__(obj.name, obj) 351 352 353class Function(NewCoreBase): 354 formula = None 355 code_string = None 356 xcode = None 357 value = None 358 _names = None 359 args = None 360 _TIME_ = None 361 362 def __init__(self, name): 363 self.setName(name) 364 self.args = [] 365 366 def __call__(self, *args): 367 for ar in range(len(args)): 368 self.__setattr__(self.args[ar], args[ar]) 369 exec(self.xcode) 370 return self.value 371 372 def setArg(self, var, value=None): 373 self.__setattr__(var, value) 374 self.args.append(var) 375 376 def addFormula(self, formula): 377 # print('PC2 377:', formula) 378 formula = formula.replace('self.', '') 379 self.formula = formula 380 InfixParser.setNameStr('self.', '') 381 InfixParser.SymbolReplacements = {'_TIME_': '_TIME_()'} 382 InfixParser.parse(formula) 383 self._names = InfixParser.names 384 self.code_string = 'self.value=%s' % InfixParser.output 385 self.xcode = compile(self.code_string, '<string>', 'exec') 386 387 388class Reaction(NewCoreBase): 389 modifiers = None 390 substrates = None 391 products = None 392 stoichiometry = None 393 multistoich = None 394 multistoich_enabled = False 395 parameters = None 396 functions = None 397 reversible = True 398 formula = None 399 code_string = None 400 rate = None 401 xcode = None 402 _names = None 403 _functions = None 404 _TIME_ = None 405 Compartment = None 406 407 def __call__(self): 408 exec(self.xcode) 409 return self.rate 410 411 def __init__(self, name): 412 self.setName(name) 413 self.modifiers = [] 414 self.substrates = [] 415 self.products = [] 416 self.stoichiometry = {} 417 self.parameters = [] 418 self.functions = [] 419 self.multistoich = [] 420 421 def addSubstrate(self, species): 422 self.__setattr__(species.name, species) 423 self.substrates.append(species) 424 425 def addProduct(self, species): 426 self.__setattr__(species.name, species) 427 self.products.append(species) 428 429 def addModifier(self, species): 430 self.__setattr__(species.name, species) 431 self.modifiers.append(species) 432 433 def addFormula(self, formula): 434 formula = formula.replace('self.', '') 435 self.formula = formula 436 InfixParser.setNameStr('self.', '()') 437 InfixParser.parse(formula) 438 self._names = InfixParser.names 439 self._functions = InfixParser.functions 440 self.code_string = 'self.rate=%s' % InfixParser.output 441 self.xcode = compile(self.code_string, '<string>', 'exec') 442 443 def addParameter(self, par): 444 self.__setattr__(par.name, par) 445 self.parameters.append(par) 446 447 def addFunction(self, func): 448 self.__setattr__(func.name, func) 449 self.functions.append(func) 450 451 def hasProducts(self, t=type): 452 return MapList([p.name for p in self.products]) 453 454 def hasSubstrates(self): 455 return MapList([s.name for s in self.substrates]) 456 457 def hasModifiers(self): 458 return MapList([m.name for m in self.modifiers]) 459 460 def hasParameters(self): 461 return MapList([p.name for p in self.parameters]) 462 463 def hasReagents(self): 464 return MapList(self.hasSubstrates() + self.hasProducts()) 465 466 def setCompartment(self, compartment): 467 self.Compartment = compartment 468 469 def getCompartment(self): 470 return self.Compartment 471 472 def hasCompartment(self): 473 if self.Compartment != None: 474 return True 475 else: 476 return False 477 478 479class Parameter(NumberBase): 480 association = None 481 482 def __init__(self, name, value): 483 self.name = name 484 self.value = value 485 self.value_initial = value 486 self.association = [] 487 488 def setAssociation(self, reac): 489 self.association.append(reac) 490 self.__setattr__(reac.name, reac) 491 492 def isParameterOf(self): 493 return MapList([a.name for a in self.association]) 494 495 496class AssignmentRule(Parameter): 497 formula = None 498 code_string = None 499 _names = None 500 _functions = None 501 type = 'assignment' 502 _TIME_ = None 503 fixed = False # added so that assignment rules can modify fixed species 504 505 def __init__(self, name, value): 506 Parameter.__init__(self, name, value) 507 508 def __call__(self): 509 exec(self.xcode) 510 return self.value 511 512 def addFormula(self, formula): 513 formula = formula.replace('self.', '') 514 self.formula = formula 515 InfixParser.setNameStr('self.', '()') 516 InfixParser.parse(formula) 517 self.code_string = 'self.value=%s' % InfixParser.output 518 self._names = InfixParser.names 519 self._functions = InfixParser.functions 520 # print(self.code_string) 521 self.xcode = compile(self.code_string, '<string>', 'exec') 522 523 def addModelAttr(self, obj): 524 self.__setattr__(obj.name, obj) 525 526 527class RateRule(NewCoreBase): 528 formula = None 529 rate = None 530 xcode = None 531 code_string = None 532 _names = None 533 _functions = None 534 compartment = None 535 536 def __init__(self, name, formula): 537 self.name = name 538 self.addFormula(formula) 539 540 def __call__(self): 541 exec(self.xcode) 542 return self.rate 543 544 def addFormula(self, formula): 545 formula = formula.replace('self.', '') 546 self.formula = formula.replace('()', '') 547 InfixParser.setNameStr('self.', '()') 548 InfixParser.parse(self.formula) 549 self.code_string = 'self.rate=%s' % InfixParser.output 550 self._names = InfixParser.names 551 self._functions = InfixParser.functions 552 self.xcode = compile(self.code_string, 'RateRule: %s' % self.name, 'exec') 553 554 def getFormula(self): 555 return self.formula 556 557 def addModelAttr(self, obj): 558 self.__setattr__(obj.name, obj) 559 560 561class ODE(NewCoreBase): 562 sdot = None 563 value = None 564 coefficients = None 565 reactions = None 566 independent = None 567 ode_terms = None 568 formula = '' 569 formula_alt = '' 570 code_string = 'self.value=' 571 code_string_alt = 'sdot=' 572 573 def __init__(self, species, independent=True): 574 self.sdot = species 575 self.name = 'ODE_' + species.name 576 self.reactions = [] 577 self.coefficients = [] 578 self.ode_terms = [] 579 self.independent = independent 580 581 def __call__(self): 582 exec(self.code_string) 583 return self.value 584 585 def addReaction(self, reaction, coefficient): 586 self.reactions.append(reaction) 587 self.coefficients.append(coefficient) 588 if coefficient > 0.0: 589 if coefficient == 1.0: 590 term = '+self.%s() ' % (reaction.name) 591 aterm = '+(%s) ' % (reaction.code_string.replace('self.rate=', '')) 592 fterm = '+%s' % (reaction.name) 593 afterm = '+ (%s) ' % (reaction.formula) 594 else: 595 term = '+%g*self.%s() ' % (abs(coefficient), reaction.name) 596 aterm = '+%g*(%s) ' % ( 597 abs(coefficient), 598 reaction.code_string.replace('self.rate=', ''), 599 ) 600 fterm = '+%g*%s' % (abs(coefficient), reaction.name) 601 afterm = '+ %g*(%s) ' % (abs(coefficient), reaction.formula) 602 else: 603 if coefficient == -1.0: 604 term = '-self.%s() ' % (reaction.name) 605 aterm = '-(%s) ' % (reaction.code_string.replace('self.rate=', '')) 606 fterm = '-%s' % (reaction.name) 607 afterm = '- (%s) ' % (reaction.formula) 608 else: 609 term = '-%g*self.%s() ' % (abs(coefficient), reaction.name) 610 aterm = '-%g*(%s) ' % ( 611 abs(coefficient), 612 reaction.code_string.replace('self.rate=', ''), 613 ) 614 fterm = '-%g*%s' % (abs(coefficient), reaction.name) 615 afterm = '- %g*(%s) ' % (abs(coefficient), reaction.formula) 616 self.ode_terms.append(term) 617 self.code_string += term 618 self.code_string_alt += aterm 619 self.formula += fterm 620 self.formula_alt += afterm 621 self.__setattr__(reaction.name, reaction) 622 623 def hasReactions(self): 624 return MapList([r.name for r in self.reactions]) 625 626 def getFormula(self): 627 return self.code_string 628 629 def getGlobalFormula(self): 630 return self.code_string_alt 631 632 633class StructMatrix(NewCoreBase): 634 """ 635 This class is specifically designed to store structural matrix information 636 give it an array and row/col index permutations it can generate its own 637 row/col labels given the label src. 638 """ 639 640 array = None 641 ridx = None 642 cidx = None 643 row = None 644 col = None 645 646 def __init__(self, array, ridx, cidx, row=None, col=None): 647 """ 648 Instantiate with array and matching row/col index arrays, optional label arrays 649 """ 650 self.array = array 651 self.ridx = ridx 652 self.cidx = cidx 653 self.row = row 654 self.col = col 655 self.shape = array.shape 656 657 def __call__(self): 658 return self.array 659 660 def getRowsByIdx(self, *args): 661 """Return the rows referenced by index (1,3,5)""" 662 return self.array.take(args, axis=0) 663 664 def getColsByIdx(self, *args): 665 """Return the columns referenced by index (1,3,5)""" 666 return self.array.take(args, axis=1) 667 668 def setRow(self, src): 669 """ 670 Assuming that the row index array is a permutation (full/subset) 671 of a source label array by supplying that source to setRow it 672 maps the row labels to ridx and creates self.row (row label list) 673 """ 674 self.row = [src[r] for r in self.ridx] 675 676 def setCol(self, src): 677 """ 678 Assuming that the col index array is a permutation (full/subset) 679 of a source label array by supplying that src to setCol 680 maps the row labels to cidx and creates self.col (col label list) 681 """ 682 self.col = [src[c] for c in self.cidx] 683 684 def getRowsByName(self, *args): 685 """Return the rows referenced by label ('s','x','d')""" 686 assert self.row != None, "\nI need row labels" 687 try: 688 return self.array.take([self.row.index(l) for l in args], axis=0) 689 except Exception as ex: 690 print(ex) 691 print("\nValid row labels are: %s" % self.row) 692 return None 693 694 def getColsByName(self, *args): 695 """Return the columns referenced by label ('s','x','d')""" 696 assert self.col != None, "\nI need column labels" 697 try: 698 return self.array.take([self.col.index(l) for l in args], axis=1) 699 except Exception as ex: 700 print(ex) 701 print("Valid column labels are: %s" % self.col) 702 return None 703 704 def getLabels(self, axis='all'): 705 """Return the matrix labels ([rows],[cols]) where axis='row'/'col'/'all'""" 706 if axis == 'row': 707 return self.row 708 elif axis == 'col': 709 return self.col 710 else: 711 return self.row, self.col 712 713 def getIndexes(self, axis='all'): 714 """Return the matrix indexes ([rows],[cols]) where axis='row'/'col'/'all'""" 715 if axis == 'row': 716 return self.ridx 717 elif axis == 'col': 718 return self.cidx 719 else: 720 return self.ridx, self.cidx 721 722 def getByIdx(self, row, col): 723 assert row in self.ridx, '\n%s is an invalid index' % row 724 assert col in self.cidx, '\n%s is an invalid index' % col 725 return self.array[row, col] 726 727 def getByName(self, row, col): 728 assert row in self.row, '\n%s is an invalid name' % row 729 assert col in self.col, '\n%s is an invalid name' % col 730 return self.array[self.row.index(row), self.col.index(col)] 731 732 def setByIdx(self, row, col, val): 733 assert row in self.ridx, '\n%s is an invalid index' % row 734 assert col in self.cidx, '\n%s is an invalid index' % col 735 self.array[row, col] = val 736 737 def setByName(self, row, col, val): 738 assert row in self.row, '\n%s is an invalid name' % row 739 assert col in self.col, '\n%s is an invalid name' % col 740 self.array[self.row.index(row), self.col.index(col)] = val 741 742 def shape(self): 743 return self.array.shape 744 745 746class EventAssignment(NumberBase): 747 variable = None 748 _names = None 749 formula = None 750 code_string = None 751 xcode = None 752 753 def __call__(self): 754 self.variable.value = self.value 755 if self.__DEBUG__: 756 print('\tAssigning %s = %s' % (self.variable.name, self.value)) 757 return True 758 759 def __init__(self, name='None'): 760 self.setName(name) 761 762 def setVariable(self, var): 763 self.variable = var 764 765 def setFormula(self, formula): 766 self.formula = formula 767 InfixParser.setNameStr('self.', '()') 768 ## InfixParser.SymbolReplacements = {'_TIME_':'_TIME_()'} 769 InfixParser.parse(formula) 770 self._names = InfixParser.names 771 self.code_string = 'self.value=%s' % InfixParser.output 772 self.xcode = compile(self.code_string, '<string>', 'exec') 773 if self.__DEBUG__: 774 '\t', self.name, self.code_string 775 776 def evaluateAssignment(self): 777 exec(self.xcode) 778 779 780class Event(NewCoreBase): 781 trigger = None 782 delay = 0.0 783 784 formula = None 785 code_string = None 786 xcode = None 787 788 state0 = False 789 state = False 790 791 assignments = None 792 _TIME_ = None 793 _ASS_TIME_ = 0.0 794 _need_action = False 795 _names = None 796 _time_symbol = None 797 798 def __init__(self, name): 799 self.setName(name) 800 self.assignments = [] 801 802 def __call__(self, time): 803 self._TIME_.set(time) 804 exec(self.xcode) 805 if self.state0 and not self.state: 806 self.state0 = self.state 807 if not self.state0 and self.state: 808 for ass in self.assignments: 809 ass.evaluateAssignment() 810 self.state0 = self.state 811 self._need_action = True 812 self._ASS_TIME_ = self._TIME_() + self.delay 813 if self.__DEBUG__: 814 print('event %s is evaluating at %s' % (self.name, time)) 815 if self._need_action and self._TIME_() >= self._ASS_TIME_: 816 for ass in self.assignments: 817 ass() 818 if self.__DEBUG__: 819 print( 820 'event %s is assigning at %s (delay=%s)' 821 % (self.name, time, self.delay) 822 ) 823 self._need_action = False 824 825 def setTrigger(self, formula, delay=0.0): 826 self.formula = formula 827 self.delay = delay 828 InfixParser.setNameStr('self.', '()') 829 ## print self._time_symbol 830 if self._time_symbol != None: 831 InfixParser.SymbolReplacements = {self._time_symbol: '_TIME_'} 832 ## self.formula = formula.replace(self._time_symbol, '_TIME_') 833 InfixParser.parse(formula) 834 self._names = InfixParser.names 835 self.code_string = 'self.state=%s' % InfixParser.output 836 if self._time_symbol != None: 837 InfixParser.setNameStr('', '') 838 InfixParser.SymbolReplacements = {self._time_symbol: '_TIME_'} 839 InfixParser.parse(formula) 840 self.formula = InfixParser.output 841 self.xcode = compile(self.code_string, '<string>', 'exec') 842 if self.__DEBUG__: 843 self.name, self.code_string 844 845 def setTriggerAttributes(self, core): 846 # TODO: experimental 847 for n in self._names: 848 self.__setattr__(n, core.__getattribute__(n)) 849 850 def setAssignment(self, var, formula): 851 ass = EventAssignment(var.name) 852 ass.setVariable(var) 853 ass.setFormula(formula) 854 self.assignments.append(ass) 855 self.__setattr__('_' + var.name, ass) 856 857 858class PieceWise(NewCoreBase): 859 """ 860 Generic piecewise class written by me! 861 862 - *args* a dictionary of piecewise information generated by the InfixParser 863 """ 864 865 name = None 866 value = None 867 formula = None 868 code_string = None 869 xcode = None 870 _names = None 871 _TIME_ = None 872 873 def __init__(self, pwd): 874 pwd = pwd.copy() 875 if pwd['other'] != None: 876 other = 'self.value = %s' % pwd.pop('other') 877 else: 878 other = 'pass' 879 pwd.pop('other') 880 InfixParser.setNameStr('self.', '') 881 InfixParser.SymbolReplacements = {'_TIME_': '_TIME_()'} 882 self._names = [] 883 if len(list(pwd.keys())) == 1: 884 formula = pwd[0][0] 885 InfixParser.parse(formula) 886 for n in InfixParser.names: 887 if n not in self._names and n != '_TIME_()': 888 self._names.append(n) 889 formula = InfixParser.output 890 self.code_string = 'if %s:\n self.value = %s\nelse:\n %s' % ( 891 formula, 892 pwd[0][1], 893 other, 894 ) 895 self.formula = self.code_string.replace('self.', '') 896 else: 897 formula = pwd[0][0] 898 InfixParser.parse(formula) 899 for n in InfixParser.names: 900 if n not in self._names and n != '_TIME_()': 901 self._names.append(n) 902 903 formula = InfixParser.output 904 self.code_string = 'if %s:\n self.value = %s\n' % (formula, pwd[0][1]) 905 pwd.pop(0) 906 for p in pwd: 907 formula = pwd[p][0] 908 InfixParser.SymbolReplacements = {'_TIME_': '_TIME_()'} 909 InfixParser.parse(formula) 910 for n in InfixParser.names: 911 if n not in self._names and n != '_TIME_()': 912 self._names.append(n) 913 formula = InfixParser.output 914 self.code_string += 'elif %s:\n self.value = %s\n' % ( 915 formula, 916 pwd[p][1], 917 ) 918 self.code_string += 'else:\n %s' % other 919 self.formula = self.code_string.replace('self.', '') 920 self.xcode = compile(self.code_string, 'PieceWise', 'exec') 921 922 def __call__(self): 923 exec(self.xcode) 924 return self.value 925 926 927class Time(object): 928 value = None 929 name = '__TIME__' 930 931 def __init__(self, t=0): 932 self.value = t 933 934 def __call__(self): 935 return self.value 936 937 def set(self, t): 938 self.value = t 939 940 941## def delay(*args): 942## print 'delay() ignored' 943## return 1.0 944 945 946class NewCore(NewCoreBase): 947 __nDict__ = None 948 reactions = None 949 species = None 950 species_variable = None 951 __model__ = None 952 __InitDict__ = None 953 __not_inited__ = None 954 global_parameters = None 955 __parameter_store__ = None 956 forcing_functions = None 957 __rules__ = None 958 __events__ = None 959 # new 960 __compartments__ = None 961 compartments = None 962 rate_rules = None 963 description = "Pysces Core2" 964 __uDict__ = None 965 stoichiometric_matrix = None 966 struct = None 967 ODEs = None 968 functions = None 969 _TIME_ = None 970 events = None 971 __sDict__ = None 972 __KeyWords__ = None 973 __piecewises__ = None 974 piecewise_functions = None 975 netStoich = None 976 977 def __init__(self, model, iValues=True, netStoich=True): 978 # setup core dictionaries 979 self.__nDict__ = model.__nDict__ 980 self.__sDict__ = model.__sDict__ 981 self.__KeyWords__ = model.__KeyWords__ 982 if self.__KeyWords__['Modelname'] != None: 983 self.setName(self.__KeyWords__['Modelname']) 984 else: 985 self.setName('PySCeSModel') 986 if self.__KeyWords__['Description'] != None: 987 self.setDescription(self.__KeyWords__['Description']) 988 else: 989 self.setDescription('PySCeSModel') 990 991 self.__model__ = model 992 self.__InitDict__ = model.__InitDict__ 993 if not iValues: 994 if self.__DEBUG__: 995 print(self.__InitDict__) 996 for k in list(self.__InitDict__.keys()): 997 self.__InitDict__[k] = getattr(self.__model__, k) 998 for c in model.__compartments__: 999 model.__compartments__[c]['size'] = getattr(self.__model__, c) 1000 self.netStoich = netStoich 1001 1002 self.global_parameters = [] 1003 self.__parameter_store__ = [] 1004 self.__not_inited__ = [] 1005 self.forcing_functions = [] 1006 self.__rules__ = model.__rules__ 1007 self.__uDict__ = model.__uDict__ 1008 self.__piecewises__ = model.__piecewises__ 1009 InfixParser.__pwcntr__ = 0 1010 1011 # start building objects 1012 self.__compartments__ = model.__compartments__ 1013 self.addCompartments() 1014 self._TIME_ = Time() 1015 self.addPieceWiseFunctions() # this adds any piecewise functions 1016 self.addSpecies() 1017 1018 # the order is important from here as eg functions can occur in rate equations 1019 try: 1020 self.__functions__ = model.__functions__ 1021 except: 1022 self.__functions__ = {} 1023 if self.__DEBUG__: 1024 print('No functions') 1025 self.functions = [] 1026 self.addFunctions() 1027 self.addReactions() 1028 self.generateMappings() 1029 self.setAssignmentRules() 1030 self.setRateRules() 1031 1032 # add event support 1033 self.__events__ = self.__model__.__eDict__ 1034 self.events = [] 1035 self.addEvents() 1036 self.addPieceWiseFunctions(update=True) # this updates their attributes 1037 1038 ## # get rid of _TIME_ in not intited 1039 ## if '_TIME_' in self.__not_inited__: 1040 ## self.__not_inited__.pop(self.__not_inited__.index('_TIME_')) 1041 if len(self.__not_inited__) > 0: 1042 print( 1043 "\nWARNING: Uninitialised parameters: %s (value set to zero)" 1044 % self.__not_inited__ 1045 ) 1046 1047 def __cleanString__(self, s): 1048 s = s.lstrip() 1049 s = s.rstrip() 1050 return s 1051 1052 ## def parseForcingFunctions(self): 1053 ## self.__rules__ = {} 1054 ## try: 1055 ## ff = self.__function_forced_str__.split('\n') 1056 ## for f in ff: 1057 ## if f != '': 1058 ## f = f.split('=') 1059 ## f[0] = f[0].replace('self.','') 1060 ## f[1] = f[1].replace('self.','') 1061 ## self.__rules__.setdefault(self.__cleanString__(f[0]), self.__cleanString__(f[1])) 1062 ## except Exception, ex: 1063 ## pass 1064 # print 'No forcing functions (%s).' % ex 1065 1066 def setDescription(self, txt): 1067 self.description = str(txt) 1068 1069 def getDescription(self): 1070 return str(self.description) 1071 1072 def setGlobalUnits(self, **kwargs): 1073 for un in list(kwargs.keys()): 1074 self.__uDict__[un] = (kwargs[un][0], kwargs[un][1]) 1075 if self.__DEBUG__: 1076 print( 1077 "Modified \"%s\" to be %i*%s*10**%i" 1078 % (un, kwargs[un][0], un, kwargs[un][1]) 1079 ) 1080 1081 def getGlobalUnits(self): 1082 return self.__uDict__ 1083 1084 def addPieceWiseFunctions(self, update=False): 1085 if not update: 1086 self.piecewise_functions = [] 1087 for pw in list(self.__piecewises__.keys()): 1088 if self.__DEBUG__: 1089 print('Info: adding piecewise function:%s' % pw) 1090 P = PieceWise(self.__piecewises__[pw]) 1091 P.setName(pw) 1092 P.__setattr__('_TIME_', self.__getattribute__('_TIME_')) 1093 self.piecewise_functions.append(P) 1094 self.__setattr__(pw, P) 1095 else: 1096 for pw in self.piecewise_functions: 1097 for a in pw._names: 1098 pw.__setattr__(a, self.__getattribute__(a)) 1099 1100 def addOneCompartment(self, name, size, dimensions, compartment=None, area=None): 1101 C = Compartment(name, compartment) 1102 C.setSize(size, dimensions) 1103 ## C.setArea(area) 1104 self.compartments.append(C) 1105 self.__setattr__(name, C) 1106 1107 def addCompartments(self): 1108 self.compartments = [] 1109 for C in self.__compartments__: 1110 c2 = self.__compartments__[C] 1111 if self.__DEBUG__: 1112 print('Adding compartment %s' % c2['name']) 1113 self.addOneCompartment( 1114 c2['name'], 1115 c2['size'], 1116 c2['dimensions'], 1117 compartment=c2['compartment'], 1118 area=None, 1119 ) 1120 1121 def addOneSpecies( 1122 self, species, value, fix=False, comp=None, amount=False, fullName=None 1123 ): 1124 s = Species(species, value) 1125 ## if comp != None: 1126 s.setCompartment(comp) 1127 s.setAmount(amount) 1128 s.setAnnotation('sbml_name', fullName) 1129 if fix: 1130 s.fixed = True 1131 self.__setattr__(species, s) 1132 self.species.append(s) 1133 if not fix: 1134 self.species_variable.append(s) 1135 if comp != None: 1136 comp.addSpecies(s) 1137 1138 def addSpecies(self): 1139 self.species = [] 1140 self.species_variable = [] 1141 for s in self.__sDict__: 1142 ## print s 1143 ## print self.__sDict__[s] 1144 name = self.__sDict__[s]['name'] 1145 if s in self.__InitDict__: 1146 val = self.__InitDict__[s] 1147 else: 1148 val = 0.0 1149 ## print val 1150 fix = self.__sDict__[s]['fixed'] 1151 if self.__sDict__[s]['compartment'] != None: 1152 comp = self.__getattribute__(self.__sDict__[s]['compartment']) 1153 else: 1154 comp = None 1155 amount = self.__sDict__[s]['isamount'] 1156 fullName = None 1157 if 'fullName' in self.__sDict__[s]: 1158 fullName = self.__sDict__[s]['fullName'] 1159 # print(self.__sDict__) 1160 self.addOneSpecies( 1161 name, val, fix=fix, comp=comp, amount=amount, fullName=fullName 1162 ) 1163 1164 def addOneFunction(self, name, args, formula): 1165 func = Function(name) 1166 # TODO: make better 1167 setattr(func, '_TIME_', self._TIME_) 1168 for a in args: 1169 func.setArg(a) 1170 func.addFormula(formula) 1171 self.functions.append(func) 1172 self.__setattr__(name, func) 1173 1174 def addFunctions(self): 1175 for f in list(self.__functions__.keys()): 1176 self.addOneFunction( 1177 f, self.__functions__[f]['args'], self.__functions__[f]['formula'] 1178 ) 1179 1180 def addOneReaction(self, rDict): 1181 r = Reaction(rDict['name']) 1182 if rDict['compartment'] != None: 1183 C = self.__getattribute__(rDict['compartment']) 1184 r.setCompartment(C) 1185 C.addReaction(r) 1186 fullName = None 1187 if 'fullName' in rDict: 1188 r.setAnnotation('sbml_name', rDict['fullName']) 1189 1190 # add dummy reaction kinetic law 1191 # TODO: make better 1192 setattr(r, '_TIME_', self._TIME_) 1193 if rDict['RateEq'] != None: 1194 r.addFormula(rDict['RateEq'].replace('self.', '')) 1195 else: 1196 r.addFormula('J') 1197 self.addParameter('J') 1198 if rDict['Type'] == 'Irrev': 1199 r.reversible = False 1200 # now we can add formulas that occured in the rate equation 1201 if len(r._functions) > 0: 1202 for func in r._functions: 1203 try: 1204 r.addFunction(self.__getattribute__(func)) 1205 except Exception as ex: 1206 print(ex) 1207 print('\nHave you added the function objects yet (addFunctions())') 1208 1209 # fxnames = self.hasFixedSpecies() 1210 processed_parameter = [] 1211 1212 # where parameters are defined `locally' per reaction 1213 for p in rDict['Params']: 1214 p = p.replace('self.', '') 1215 if p not in self.hasGlobalParameters() and not ( 1216 p in self.hasFixedSpecies() or p in self.__compartments__ 1217 ): 1218 if self.__DEBUG__: 1219 print("Adding parameter %s from networkdict" % p) 1220 self.addParameter(p) 1221 par = self.__getattribute__(p) 1222 par.setAssociation(r) 1223 r.addParameter(par) 1224 processed_parameter.append(p) 1225 elif not (p in self.hasFixedSpecies() or p in self.__compartments__): 1226 if self.__DEBUG__: 1227 print("Updating parameter %s from networkdict" % p) 1228 pidx = self.hasGlobalParameters().index(p) 1229 self.global_parameters[pidx].setAssociation(r) 1230 r.addParameter(self.global_parameters[pidx]) 1231 processed_parameter.append(p) 1232 1233 # print self.hasGlobalParameters() 1234 # where parameters are not `locally' defined and are extracted from Req (ie from SBML) 1235 for p in r._names: 1236 p = p.replace('self.', '') 1237 if p == '_TIME_': 1238 pass 1239 elif p in [pw.name for pw in self.piecewise_functions]: 1240 pass 1241 elif p in self.hasCompartments() and p not in processed_parameter: 1242 C = self.__getattribute__(p) 1243 C.addReaction(r) 1244 # TODO: this will work until isParameterOf is called on a compartment object 1245 r.addParameter(C) 1246 # dirty alternative 1247 # setattr(r, C.name, C) 1248 processed_parameter.append(p) 1249 elif ( 1250 p not in processed_parameter 1251 and p not in self.hasGlobalParameters() 1252 and p not in self.hasSpecies() 1253 ): 1254 if self.__DEBUG__: 1255 print("Adding parameter %s from global" % p) 1256 self.addParameter(p) 1257 par = self.__getattribute__(p) 1258 par.setAssociation(r) 1259 r.addParameter(par) 1260 processed_parameter.append(p) 1261 1262 elif p not in processed_parameter and p not in self.hasSpecies(): 1263 if self.__DEBUG__: 1264 print("Updating parameter %s from global" % p) 1265 pidx = self.hasGlobalParameters().index(p) 1266 self.global_parameters[pidx].setAssociation(r) 1267 r.addParameter(self.global_parameters[pidx]) 1268 processed_parameter.append(p) 1269 1270 self.__setattr__(rDict['name'], r) 1271 self.reactions.append(r) 1272 1273 def addParameter(self, name): 1274 if name not in self.__piecewises__: 1275 if name in self.__InitDict__: 1276 par = Parameter(name, self.__InitDict__[name]) 1277 else: 1278 par = Parameter(name, 0.0) 1279 if name not in self.__not_inited__: 1280 self.__not_inited__.append(name) 1281 self.global_parameters.append(par) 1282 self.__setattr__(name, par) 1283 1284 def addReactions(self): 1285 self.reactions = [] 1286 for r in self.__model__.reactions: 1287 self.addOneReaction(self.__nDict__[r]) 1288 non_parameters = ( 1289 self.hasGlobalParameters() + self.hasSpecies() + self.hasFixedSpecies() 1290 ) 1291 for k in list(self.__InitDict__.keys()): 1292 if k not in non_parameters: 1293 if self.__DEBUG__: 1294 print('Adding new parameter:', k) 1295 self.addParameter(k) 1296 1297 def replaceParameterWithRule(self, ar): 1298 par = self.__getattribute__(ar.name) 1299 for r in par.association: 1300 ar.setAssociation(r) 1301 setattr(r, ar.name, ar) 1302 r.parameters[r.hasParameters().index(ar.name)] = ar 1303 self.global_parameters[self.hasGlobalParameters().index(ar.name)] = ar 1304 self.__setattr__(ar.name, ar) 1305 1306 def replaceFixedSpeciesWithRule(self, ar): 1307 fs = self.__getattribute__(ar.name) 1308 ar.fixed = fs.fixed 1309 for r in fs.subs: 1310 ar.setSubstrate(r) 1311 setattr(r, ar.name, ar) 1312 r.substrates[r.hasSubstrates().index(ar.name)] = ar 1313 for r in fs.prods: 1314 ar.setProduct(r) 1315 setattr(r, ar.name, ar) 1316 r.products[r.hasProducts().index(ar.name)] = ar 1317 for r in fs.mods: 1318 ar.setModifier(r) 1319 setattr(r, ar.name, ar) 1320 r.modifiers[r.hasModifiers().index(ar.name)] = ar 1321 self.species[self.hasSpecies().index(ar.name)] = ar 1322 self.__setattr__(ar.name, ar) 1323 1324 def replaceSpeciesWithRule(self, ar): 1325 fs = self.__getattribute__(ar.name) 1326 for r in fs.subs: 1327 ar.setSubstrate(r) 1328 setattr(r, ar.name, ar) 1329 r.substrates[r.hasSubstrates().index(ar.name)] = ar 1330 for r in fs.prods: 1331 ar.setProduct(r) 1332 setattr(r, ar.name, ar) 1333 r.products[r.hasProducts().index(ar.name)] = ar 1334 for r in fs.mods: 1335 ar.setModifier(r) 1336 setattr(r, ar.name, ar) 1337 r.modifiers[r.hasModifiers().index(ar.name)] = ar 1338 self.species[self.hasSpecies().index(ar.name)] = ar 1339 self.species_variable[self.hasVariableSpecies().index(ar.name)] = ar 1340 self.__setattr__(ar.name, ar) 1341 1342 def setAssignmentRules(self): 1343 aps = [ 1344 self.__rules__[ar]['name'] 1345 for ar in self.__rules__ 1346 if self.__rules__[ar]['type'] == 'assignment' 1347 ] 1348 ## for p in self.global_parameters + [self.get(fs) for fs in self.hasFixedSpecies()]: 1349 for p in self.global_parameters + self.species: 1350 # print p.name 1351 if p.name in aps: 1352 if self.__DEBUG__: 1353 print( 1354 'Assigning: %s = %s' 1355 % (p.name, self.__rules__[p.name]['formula']) 1356 ) 1357 p2 = None 1358 # TODO: make better 1359 if p.name in self.hasGlobalParameters(): 1360 p2 = AssignmentRule(p.name, self.__InitDict__[p.name]) 1361 setattr(p2, '_TIME_', self._TIME_) 1362 self.replaceParameterWithRule(p2) 1363 elif p.name in self.hasFixedSpecies(): 1364 p2 = SpeciesAssignmentRule(p.name, self.__InitDict__[p.name]) 1365 p2.setCompartment(p.getCompartment()) 1366 setattr(p2, '_TIME_', self._TIME_) 1367 self.replaceFixedSpeciesWithRule(p2) 1368 elif p.name in self.hasVariableSpecies(): 1369 p2 = SpeciesAssignmentRule(p.name, self.__InitDict__[p.name]) 1370 p2.setCompartment(p.getCompartment()) 1371 setattr(p2, '_TIME_', self._TIME_) 1372 self.replaceSpeciesWithRule(p2) 1373 1374 assert isinstance(p2, AssignmentRule) or isinstance( 1375 p2, SpeciesAssignmentRule 1376 ), "\nHappy assertion error" 1377 # print type(p2) 1378 p2.addFormula(self.__rules__[p.name]['formula']) 1379 1380 # add piecewise functions to assignment rules or something like that 1381 if 'piecewise' in p2.formula: 1382 print('\n\nWe have a piecewise', p2.name, p2.formula) 1383 print(p2, self.piecewise_functions[0]._names[0], self.__getattribute__(self.piecewise_functions[0]._names[0]).__call__()) 1384 for p in range(len(self.piecewise_functions)): 1385 for n in range(len(self.piecewise_functions[p]._names)): 1386 setattr(p2, self.piecewise_functions[p]._names[n], self.__getattribute__(self.piecewise_functions[p]._names[n]).__call__()) 1387 self.piecewise_functions[p].code_string = self.piecewise_functions[p].code_string.replace('self.{}'.format(self.piecewise_functions[p]._names[n]), 'self.{}()'.format(self.piecewise_functions[p]._names[n])) 1388 self.piecewise_functions[p].xcode = compile(self.piecewise_functions[p].code_string, 'PieceWise', 'exec') 1389 setattr(p2, self.piecewise_functions[p].name, self.piecewise_functions[p]) 1390 print(self.piecewise_functions[p].code_string) 1391 1392 1393 for n in p2._names + p2._functions: 1394 p2.addModelAttr(self.__getattribute__(n)) 1395 ## # setup initial values 1396 ## p2.value_initial = self.p2() 1397 if p2.name in self.__not_inited__: 1398 self.__not_inited__.pop(self.__not_inited__.index(p.name)) 1399 for p in self.global_parameters: 1400 if p.name in self.hasAssignmentRules(): 1401 # TODO assignment rules need a list of properties 1402 for ar in p._names: 1403 if ar in self.hasAssignmentRules(): 1404 setattr(p, ar, self.__getattribute__(ar)) 1405 # TODO this is where things will go wrong if fs --> ar contains nested ar's 1406 1407 def setRateRules(self): 1408 # TODO mayvbe split into two methods for now read from self.__rules__ 1409 # TODO add functions to rules 1410 ars = [ 1411 self.__rules__[ar]['name'] 1412 for ar in self.__rules__ 1413 if self.__rules__[ar]['type'] == 'rate' 1414 ] 1415 self.rate_rules = [] 1416 for rr in ars: 1417 rrobj = RateRule(self.__rules__[rr]['name'], self.__rules__[rr]['formula']) 1418 ## print 'RR:', rrobj.name, rrobj._names, rrobj._functions 1419 for symb in rrobj._names + rrobj._functions: 1420 rrobj.addModelAttr(self.__getattribute__(symb)) 1421 self.rate_rules.append(rrobj) 1422 # TODO investgiate this as it is problematic, the rate rule 1423 # is not a model property as such more an ODE property 1424 ## self.__setattr__(rrobj.name, rrobj) 1425 if self.__DEBUG__: 1426 print( 1427 'Adding RateRule %s with formula: %s' % (rrobj.name, rrobj.formula) 1428 ) 1429 1430 def addOneEvent(self, e): 1431 """Add a single event using an event dictionary """ 1432 # translate self.__events__[e] to e 1433 1434 ev = Event(e['name']) 1435 ev._time_symbol = e['tsymb'] 1436 ev.setTrigger(e['trigger'], e['delay']) 1437 # associate model attributes with event 1438 # TODO: check that this still works 1439 ev.setTriggerAttributes(self) 1440 ## for n in ev._names: 1441 ## setattr(ev, n, self.__getattribute__(n)) 1442 # for each assignment 1443 for ass in e['assignments']: 1444 ev.setAssignment(self.__getattribute__(ass), e['assignments'][ass]) 1445 assref = getattr(ev, '_' + ass) # don\t like this at all :-( 1446 # associate model attributes with assignment 1447 for n in assref._names: 1448 setattr(assref, n, self.__getattribute__(n)) 1449 self.events.append(ev) 1450 self.__setattr__(ev.name, ev) 1451 setattr(ev, '_TIME_', self._TIME_) 1452 1453 def addEvents(self): 1454 # TODO: check that you can change the trigger on the fly (might need a setAttr thing in event obj) 1455 self.events = [] 1456 # for each event 1457 for e in self.__events__: 1458 self.addOneEvent(self.__events__[e]) 1459 1460 def generateMappings(self): 1461 ## self.netStoich = False 1462 for reac in self.reactions: 1463 if self.netStoich: 1464 for reag in self.__nDict__[reac.name]['Reagents']: 1465 if self.__nDict__[reac.name]['Reagents'][reag] < 0.0: 1466 reac.addSubstrate( 1467 self.__getattribute__(reag.replace('self.', '')) 1468 ) 1469 self.__getattribute__(reag.replace('self.', '')).setSubstrate( 1470 self.__getattribute__(reac.name) 1471 ) 1472 else: 1473 reac.addProduct( 1474 self.__getattribute__(reag.replace('self.', '')) 1475 ) 1476 self.__getattribute__(reag.replace('self.', '')).setProduct( 1477 self.__getattribute__(reac.name) 1478 ) 1479 reac.stoichiometry.setdefault( 1480 reag.replace('self.', ''), 1481 self.__nDict__[reac.name]['Reagents'][reag], 1482 ) 1483 else: 1484 for reag in self.__nDict__[reac.name]['AllReagents']: 1485 if reag[1] < 0.0: 1486 reac.addSubstrate( 1487 self.__getattribute__(reag[0].replace('self.', '')) 1488 ) 1489 self.__getattribute__( 1490 reag[0].replace('self.', '') 1491 ).setSubstrate(self.__getattribute__(reac.name)) 1492 else: 1493 reac.addProduct( 1494 self.__getattribute__(reag[0].replace('self.', '')) 1495 ) 1496 self.__getattribute__(reag[0].replace('self.', '')).setProduct( 1497 self.__getattribute__(reac.name) 1498 ) 1499 reac.multistoich.append((reag[0].replace('self.', ''), reag[1])) 1500 if reag[0].replace('self.', '') in reac.stoichiometry: 1501 reac.multistoich_enabled = True 1502 reac.stoichiometry.setdefault(reag[0].replace('self.', ''), reag[1]) 1503 for mod in self.__nDict__[reac.name]['Modifiers']: 1504 reac.addModifier(self.__getattribute__(mod.replace('self.', ''))) 1505 self.__getattribute__(mod.replace('self.', '')).setModifier( 1506 self.__getattribute__(reac.name) 1507 ) 1508 ## print 'I AM LEGEND' 1509 ## print reac.stoichiometry 1510 ## print reac.multistoich 1511 ## print 'reac.multistoich_enabled', reac.multistoich_enabled 1512 ## print self.__nDict__[reac.name]['Reagents'] 1513 ## print self.__nDict__[reac.name]['AllReagents'] 1514 1515 def setStoichiometricMatrix(self): 1516 vspec = self.hasVariableSpecies() 1517 react = self.hasReactions() 1518 nm = numpy.zeros((len(vspec), len(react)), 'd') 1519 for sp in vspec: 1520 for r in self.get(sp).isReagentOf(): 1521 nm[vspec.index(sp)][react.index(r)] = self.get(r).stoichiometry[sp] 1522 # this is if absolute stoichiometry value is used 1523 ## for r in self.get(sp).isSubstrateOf(): 1524 ## nm[vspec.index(sp)][react.index(r)] = abs(self.get(r).stoichiometry[sp]) 1525 ## for r in self.get(sp).isProductOf(): 1526 ## nm[vspec.index(sp)][react.index(r)] = -abs(self.get(r).stoichiometry[sp]) 1527 self.stoichiometric_matrix = StructMatrix( 1528 nm, list(range(len(vspec))), list(range(len(react))) 1529 ) 1530 self.stoichiometric_matrix.setRow(vspec) 1531 self.stoichiometric_matrix.setCol(react) 1532 1533 def addODEs(self): 1534 self.ODEs = [] 1535 for varspec in self.stoichiometric_matrix.row: 1536 if self.struct != None: 1537 if varspec not in self.struct.Nr.row: 1538 if self.__DEBUG__: 1539 print('Creating dependent ODE_%s' % varspec) 1540 ode = ODE(self.get(varspec), independent=False) 1541 else: 1542 if self.__DEBUG__: 1543 print('Creating independent ODE_%s' % varspec) 1544 ode = ODE(self.get(varspec), independent=True) 1545 else: 1546 if self.__DEBUG__: 1547 print( 1548 'Creating independent* ODE_%s (*assumed - no structural information available)' 1549 % varspec 1550 ) 1551 ode = ODE(self.get(varspec), independent=True) 1552 mrow = self.stoichiometric_matrix.getRowsByName(varspec) 1553 for e in range(len(mrow[0])): 1554 if mrow[0, e] != 0.0: 1555 print( 1556 'Adding term: %s*%s' 1557 % (mrow[0, e], self.stoichiometric_matrix.col[e]) 1558 ) 1559 ode.addReaction( 1560 self.get(self.stoichiometric_matrix.col[e]), mrow[0, e] 1561 ) 1562 self.__setattr__(ode.name, ode) 1563 self.ODEs.append(ode) 1564 self.__setattr__( 1565 'xcode_' + ode.name, compile(ode.getGlobalFormula(), '<string>', 'exec') 1566 ) 1567 1568 def hasODEs(self): 1569 return MapList([o.name for o in self.ODEs]) 1570 1571 def evalODEs(self, odes): 1572 return [v() for v in odes] 1573 1574 def evalXcode(self, ode): 1575 exec(self.__getattribute__('xcode_' + ode.name)) 1576 return sdot 1577 1578 def hasFunctions(self): 1579 return MapList([f.name for f in self.functions]) 1580 1581 def hasReactions(self): 1582 return MapList([r.name for r in self.reactions]) 1583 1584 def hasSpecies(self): 1585 return MapList([s.name for s in self.species]) 1586 1587 def hasFixedSpecies(self): 1588 return MapList([s.name for s in self.species if s.fixed]) 1589 1590 def hasVariableSpecies(self): 1591 return MapList([s.name for s in self.species if not s.fixed]) 1592 1593 def findReactionsThatIncludeAllSpecifiedReagents(self, *args): 1594 assert len(args) > 1, '\nNeed two or more species for this one!' 1595 setlist = [self.__getattribute__(s).isReagentOf().asSet() for s in args] 1596 isect = setlist[0] 1597 for s in setlist: 1598 isect.intersection_update(s) 1599 return MapList(isect) 1600 1601 def hasGlobalParameters(self): 1602 return MapList(p.name for p in self.global_parameters) 1603 1604 def hasAssignmentRules(self): 1605 return MapList( 1606 [ 1607 ar.name 1608 for ar in self.global_parameters + self.species 1609 if hasattr(ar, 'type') == 'assignment' 1610 ] 1611 ) 1612 1613 #def hasAssignmentRules(self): 1614 #return MapList( 1615 #[ 1616 #ar.name 1617 #for ar in self.global_parameters + self.species 1618 #if hasattr(ar, 'type') == 'rate' 1619 #] 1620 #) 1621 1622 def hasEvents(self): 1623 return MapList(e.name for e in self.events) 1624 1625 def hasCompartments(self): 1626 return MapList(c.name for c in self.compartments) 1627 1628 1629## if __psyco_active__: 1630## psyco.bind(NewCoreBase) 1631## psyco.bind(NumberBase) 1632## psyco.bind(Species) 1633## psyco.bind(Parameter) 1634## psyco.bind(AssignmentRule) 1635## psyco.bind(Reaction) 1636## psyco.bind(ODE) 1637## psyco.bind(NewCore) 1638