SloppyCell.ReactionNetworks.Network

1 # This makes integer like one would expect mathematically, e.g. 1/2 = .5 2 # rather than 0. This is obviously safer when loading other folks' models. 3 # It does, however, cost about 10% performance on PC12. 4 # This is expected to become standard behavior around python 3.0 5 from __future__ import division 6 7 import copy 8 import sets 9 import types 10 import time 11 import os 12 import sys 13 import operator 14 15 import logging 16 logger = logging.getLogger('ReactionNetworks.Network_mod') 17 18 import scipy 19 import math 20 21 import SloppyCell 22 import SloppyCell.Utility as Utility 23 import SloppyCell.KeyedList_mod 24 KeyedList = SloppyCell.KeyedList_mod.KeyedList 25 26 import SloppyCell.ExprManip as ExprManip 27 # We load a dictionary of previously-taken derivatives for efficiency 28 ExprManip.load_derivs(os.path.join(SloppyCell._TEMP_DIR, 'diff.pickle')) 29 # This will save the diffs dictionary upon exit from the python interpreter, 30 # but only on the master node. 31 if SloppyCell.my_rank == 0: 32 import atexit 33 atexit.register(ExprManip.save_derivs, os.path.join(SloppyCell._TEMP_DIR, 34 'diff.pickle')) 35 36 import Reactions 37 from Components import * 38 import Dynamics 39 import Trajectory_mod 40 41 _double_max_ = scipy.finfo(scipy.float_).max 42 _double_tiny_ = scipy.finfo(scipy.float_).tiny 43

44 -class Network:

45 # These are the methods that should be called to generate all the dynamic 46 # functions. Each method should store the resulting function body in 47 # self._dynamic_funcs_python or self._dynamic_funcs_c. C functions should 48 # also have prototypes stored in self._prototypes_c. 49 _dynamic_structure_methods = ['_make_res_function', 50 '_make_alg_deriv_func', 51 '_make_alg_res_func', 52 '_make_dres_dc_function', 53 '_make_dres_dcdot_function', 54 '_make_ddaskr_jac', 55 '_make_dres_dsinglep', 56 '_make_sens_rhs', 57 '_make_log_funcs', 58 '_make_integrate_stochastic_tidbit' 59 ] 60 _dynamic_event_methods = ['_make_root_func'] 61 62 # These are-predefined functions we want in our working namespace 63 _common_namespace = {'log': scipy.log, 64 'log10': scipy.log10, 65 'exp': math.exp, 66 'cos': math.cos, 67 'sin': math.sin, 68 'tan': math.tan, 69 'acos': math.acos, 70 'asin': math.asin, 71 'atan': math.atan, 72 'cosh': math.cosh, 73 'sinh': math.sinh, 74 'tanh': math.tanh, 75 # These don't have C support, so I've dropped them for 76 # now. 77 #'arccosh': scipy.arccosh, 78 #'arcsinh': scipy.arcsinh, 79 #'arctanh': scipy.arctanh, 80 'pow': math.pow, 81 'sqrt': math.sqrt, 82 'exponentiale': math.e, 83 'pi': math.pi, 84 'scipy': scipy, 85 'operator': operator, 86 } 87 # These are functions we need to create but that should be used commonly 88 _standard_func_defs = [('root', ['n', 'x'], 'x**(1./n)'), 89 ('cot', ['x'], '1./tan(x)'), 90 ('arccot', ['x'], 'atan(1/x)'), 91 ('coth', ['x'], '1./tanh(x)'), 92 #('arccoth', ['x'], 'arctanh(1./x)'), 93 ('csc', ['x'], '1./sin(x)'), 94 ('arccsc', ['x'], 'asin(1./x)'), 95 ('csch', ['x'], '1./sinh(x)'), 96 #('arccsch', ['x'], 'arcsinh(1./x)'), 97 ('sec', ['x'], '1./cos(x)'), 98 ('arcsec', ['x'], 'acos(1./x)'), 99 ('sech', ['x'], '1./cosh(x)'), 100 #('arcsech', ['x'], 'arccosh(1./x)'), 101 ] 102 # These define all the functions needed to deal with triggers. 103 _logical_comp_func_defs = [('gt', ['x', 'y'], 'x > y'), 104 ('geq', ['x', 'y'], 'x >= y'), 105 ('lt', ['x', 'y'], 'x < y'), 106 ('leq', ['x', 'y'], 'x <= y'), 107 ('eq', ['x', 'y'], 'y == x'), 108 ('not_func', ['x'], 'not x'), 109 ('and_func', '*', 'x and y'), 110 ('or_func', '*', 'x or y'), 111 ] 112 # Note that and_func and or_func are special-cased below, because they can 113 # take a variable number of arguments in SBML. This has to be done because 114 # the comp_func_defs list is used two ways. 1) It is used to substitute 115 # within expressions derived from libsbml. 2) It is used to generate lambda 116 # expressions for the functions. 117 118 # Add our common function definitions to a func_strs list. 119 _common_func_strs = KeyedList() 120 for id, vars, math in _standard_func_defs: 121 var_str = ','.join(vars) 122 func = 'lambda %s: %s' % (var_str, math) 123 _common_func_strs.set(id, func) 124 # These are all the partial derivatives 125 for ii, wrt in enumerate(vars): 126 deriv_id = '%s_%i' % (id, ii) 127 diff_math = ExprManip.diff_expr(math, wrt) 128 func = 'lambda %s: %s' % (var_str, diff_math) 129 _common_func_strs.set(deriv_id, func) 130 131 # Now add the C versions. 132 _common_func_strs_c = KeyedList() 133 _common_prototypes_c = KeyedList() 134 for id, vars, math in _standard_func_defs: 135 vars_c = ['double %s' % var for var in vars] 136 var_str_c = ','.join(vars_c) 137 c_body = [] 138 c_body.append('double %s(%s){' % (id, var_str_c)) 139 c_body.append('return %s;' % ExprManip.make_c_compatible(math)) 140 c_body.append('}') 141 c_body = os.linesep.join(c_body) 142 _common_prototypes_c.set(id, 'double %s(%s);' % (id, var_str_c)) 143 _common_func_strs_c.set(id, c_body) 144 for ii, wrt in enumerate(vars): 145 deriv_id = '%s_%i' % (id, ii) 146 diff_math = ExprManip.diff_expr(math, wrt) 147 c_body = [] 148 c_body.append('double %s(%s){' % (deriv_id, var_str_c)) 149 c_body.append('return %s;' % ExprManip.make_c_compatible(diff_math)) 150 c_body.append('}') 151 c_body = os.linesep.join(c_body) 152 _common_prototypes_c.set(deriv_id,'double %s(%s);' % (deriv_id, 153 var_str_c)) 154 _common_func_strs_c.set(deriv_id, c_body) 155 156 # Also do the logical functions. These don't have derivatives. 157 for id, vars, math in _logical_comp_func_defs: 158 # and_func and or_func are special-cased, because the SBML versions can 159 # take multiple arguments. 160 if id == 'and_func': 161 func = 'lambda *x: reduce(operator.and_, x)' 162 elif id == 'or_func': 163 func = 'lambda *x: reduce(operator.or_, x)' 164 else: 165 var_str = ','.join(vars) 166 func = 'lambda %s: %s' % (var_str, math) 167 _common_func_strs.set(id, func) 168 169 # These are the strings we eval to create the extra functions our 170 # Network defines. We'll evaluate them all to start with, to get the 171 # common ones in there. 172 for func_id, func_str in _common_func_strs.items(): 173 _common_namespace[func_id] = eval(func_str, _common_namespace, {}) 174

175 - def __init__(self, id, name=''):

176 self.id, self.name = id, name 177 178 self.functionDefinitions = KeyedList() 179 self.reactions = KeyedList() 180 self.assignmentRules = KeyedList() 181 self.rateRules = KeyedList() 182 self.algebraicRules = KeyedList() 183 self.events = KeyedList() 184 self.constraints = KeyedList() 185 186 # Variables is primary storage for all compartments, species, and 187 # parameters. We build up 'cross reference' lists for convenience. 188 self.variables = KeyedList() 189 190 # The following are all 'cross reference' lists which are intended 191 # only to hold references to objects in self.variables 192 # (_makeCrossReferences will fill them in) 193 self.assignedVars = KeyedList() 194 self.constantVars = KeyedList() 195 self.optimizableVars = KeyedList() 196 self.dynamicVars = KeyedList() 197 self.algebraicVars = KeyedList() 198 self.compartments = KeyedList() 199 self.parameters = KeyedList() 200 self.species = KeyedList() 201 202 # All expressions are evaluated, and functions exec'd in self.namespace. 203 # (A dictionary mapping names to objects they represent) 204 self.namespace = copy.copy(self._common_namespace) 205 self._func_strs = copy.copy(self._common_func_strs) 206 207 # Option to integrate with log concentrations (to avoid negative 208 # concentrations) 209 self.integrateWithLogs = False 210 211 self.len_root_func = 0 212 213 # These dictionaries are for storing the function bodies of our dynamic 214 # functions. 215 self._dynamic_funcs_python = KeyedList() 216 self._prototypes_c = self._common_prototypes_c.copy() 217 self._dynamic_funcs_c = KeyedList() 218 self._func_defs_c = KeyedList() 219 220 self.deriv_funcs_enabled = True 221 self.compiled = False

222 223 add_int_times, add_tail_times = True, True

224 - def full_speed(cls):

225 """ 226 Do not add timepoints to returned trajectories. 227 """ 228 cls.add_int_times, cls.add_tail_times = False, False

229 full_speed = classmethod(full_speed)

230 - def fill_traj(cls):

231 """ 232 Allow the integrator to automatically add timepoints to the interior of 233 a trajectory. 234 235 This is slower, but may capture the dynamics better. 236 """ 237 cls.add_int_times, cls.add_tail_times = True, False

238 fill_traj = classmethod(fill_traj)

239 - def pretty_plotting(cls):

240 """ 241 Add timepoints to the interior of a trajectory, plus 5% past the last 242 timepoint requested. 243 244 This produces the prettiest plots. 245 """ 246 cls.add_int_times, cls.add_tail_times = True, True

247 pretty_plotting = classmethod(pretty_plotting) 248 249 # 250 # Methods used to build up a network 251 #

252 - def _add_variable(self, var):

253 self._checkIdUniqueness(var.id) 254 self.variables.set(var.id, var) 255 self._makeCrossReferences()

256

257 - def add_compartment(self, id, initial_size=1.0, name='', 258 typical_value=None, 259 is_constant=True, is_optimizable=False):

260 """ 261 Add a compartment to the Network. 262 263 All species must reside within a compartment. 264 """ 265 compartment = Compartment(id, initial_size, name, typical_value, 266 is_constant, is_optimizable) 267 self._add_variable(compartment)

268

269 - def add_species(self, id, compartment, initial_conc=0, 270 name='', typical_value=None, 271 is_boundary_condition=False, is_constant=False, 272 is_optimizable=False, uniprot_ids=set()):

273 """ 274 Add a species to the Network. 275 """ 276 species = Species(id, compartment, initial_conc, name, typical_value, 277 is_boundary_condition, is_constant, is_optimizable, 278 uniprot_ids) 279 self._add_variable(species)

280

281 - def add_parameter(self, id, initial_value=1.0, name='', 282 typical_value=None, 283 is_constant=True, is_optimizable=True):

284 """ 285 Add a parameter to the Network. 286 """ 287 parameter = Parameter(id, initial_value, name, is_constant, 288 typical_value, is_optimizable) 289 self._add_variable(parameter)

290

291 - def add_event(self, id, trigger, event_assignments={}, delay=0, name='', 292 buffer=0):

293 """ 294 Add an event to the Network. 295 296 id - id for this event 297 trigger - The event firest when trigger passes from False to True. 298 Examples: To fire when time becomes greater than 5.0: 299 trigger = 'gt(time, 5.0)' 300 To fire when A becomes less than sin(B/C): 301 trigger = 'lt(A, sin(B/C))' 302 event_assignments - A dictionary or KeyedList of assignments to make 303 when the event executes. 304 Example: To set A to 4.3 and D to B/C 305 event_assignments = {'A': 4.3, 306 'D': 'B/C'} 307 delay - Optionally, assignments may take effect some time after the 308 event fires. delay may be a number or math expression 309 name - A more detailed name for the event, not restricted to the id 310 format 311 """ 312 event = Event(id, trigger, event_assignments, delay, name, 313 buffer) 314 self._checkIdUniqueness(event.id) 315 self.events.set(event.id, event)

316

317 - def add_constraint(self, id, trigger, message=None, name=''):

318 """ 319 Add a constraint to the Network. 320 321 id - id for this constraint 322 trigger - We treat constraints as events that correspond to an 323 invalid solution whenever the trigger is True. 324 Example: To have an invalid solution when species A is > 5.0: 325 trigger = 'gt('A', 5.0)' 326 name - A more detailed name for the constraint, not restricted to the id 327 format 328 """ 329 constraint = ConstraintEvent(id, trigger, message, name) 330 self._checkIdUniqueness(constraint.id) 331 self.constraints.set(constraint.id, constraint)

332 333

334 - def add_func_def(self, id, variables, math, name=''):

335 """ 336 Add a function definition to the Network. 337 338 id - id for the function definition 339 variables - The variables used in the math expression whose 340 values should be subsituted. 341 math - The math expression the function definition represents 342 name - A more extended name for the definition 343 344 Example: 345 To define f(x, y, z) = y**2 - cos(x/z) 346 net.add_func_def('my_func', ('x', 'y', 'z'), 'y**2 - cos(x/z)') 347 """ 348 func = FunctionDefinition(id, variables, math, name) 349 self._checkIdUniqueness(func.id) 350 self.functionDefinitions.set(func.id, func) 351 352 # Add the function and its partial derivatives to func_strs 353 # Also do the evaluation. 354 var_str = ','.join(variables) 355 func_str = 'lambda %s: %s' % (var_str, math) 356 self._func_strs.set(id, func_str) 357 self.namespace[id] = eval(func_str, self.namespace)

358

359 - def make_func_defs(self):

360 self.namespace = copy.copy(self._common_namespace) 361 for id, func in self.functionDefinitions.items(): 362 variables = func.variables 363 math = func.math 364 365 var_str = ','.join(variables) 366 func_str = 'lambda %s: %s' % (var_str, math) 367 self._func_strs.set(id, func_str) 368 logger.debug('Adding function to namespace: %s, %s' % (id, func_str)) 369 self.namespace[id] = eval(func_str, self.namespace) 370 371 # Add the function and its partial derivatives to the C code. 372 c_vars = ['double %s' % var for var in variables] 373 c_var_str = ','.join(c_vars) 374 self._prototypes_c.set(id,'double %s(%s);' % (id, c_var_str)) 375 c_math = ExprManip.make_c_compatible(math) 376 c_body = 'double %s(%s){return %s;}' % (id, c_var_str, c_math) 377 self._func_defs_c.set(id, c_body) 378 379 for ii, wrt in enumerate(variables): 380 # Python derivatives 381 diff_id = '%s_%i' % (id, ii) 382 diff_math = ExprManip.diff_expr(math, wrt) 383 func_str = 'lambda %s: %s' % (var_str, diff_math) 384 self._func_strs.set(diff_id, func_str) 385 self.namespace[diff_id] = eval(func_str, self.namespace) 386 387 # C derivatives 388 self._prototypes_c.set(diff_id, 389 'double %s(%s);' % (diff_id,c_var_str)) 390 c_math = ExprManip.make_c_compatible(diff_math) 391 c_body = 'double %s(%s){return %s;}' % (diff_id, c_var_str, 392 c_math) 393 self._func_defs_c.set(diff_id, c_body)

394

395 - def addReaction(self, id, *args, **kwargs):

396 # Reactions can be added by (1) passing in a string representing 397 # kinetic law, or (2) passing in a class already specifying the 398 # kinetic law. 399 # XXX: I'm a little unhappy with this because option (2) breaks the 400 # pattern that the first argument is the id 401 if isinstance(id, str): 402 rxn = apply(Reactions.Reaction, (id,) + args, kwargs) 403 else: 404 rxn = apply(id, args, kwargs) 405 406 self._checkIdUniqueness(rxn.id) 407 self.reactions.set(rxn.id, rxn)

408

409 - def add_assignment_rule(self, var_id, rhs):

410 """ 411 Add an assignment rule to the Network. 412 413 A rate rules species that <var_id> = rhs. 414 """ 415 self.set_var_constant(var_id, False) 416 self.assignmentRules.set(var_id, rhs) 417 self._makeCrossReferences() 418 # Put this assignment into effect... 419 # The time=0 is somewhat arbitrary, reflecting the fact that we usually 420 # start integrations from time=0. 421 self.updateAssignedVars(time=0)

422

423 - def add_rate_rule(self, var_id, rhs):

424 """ 425 Add a rate rule to the Network. 426 427 A rate rules species that d <var_id>/dt = rhs. 428 """ 429 self.set_var_constant(var_id, False) 430 self.rateRules.set(var_id, rhs) 431 self._makeCrossReferences()

432

433 - def add_algebraic_rule(self, rhs):

434 """ 435 Add an algebraic rule to the Network. 436 437 An algebraic rule specifies that 0 = rhs. 438 """ 439 self.algebraicRules.set(rhs, rhs) 440 self._makeCrossReferences()

441

442 - def remove_component(self, id):

443 """ 444 Remove the component with the given id from the Network. 445 446 Components that can be removed are variables, reactions, events, 447 function definitions, assignment rules, rate rules, and constraints. 448 """ 449 complists = [self.variables, self.reactions, self.functionDefinitions, 450 self.events, self.constraints, self.assignmentRules, self.rateRules, 451 self.algebraicRules] 452 for complist in complists: 453 # If the id is in a list and has a non-empty name 454 if complist.has_key(id): 455 complist.remove_by_key(id) 456 457 self._makeCrossReferences()

458

459 - def _checkIdUniqueness(self, id):

460 """ 461 Check whether a given id is already in use by this Network. 462 """ 463 if id == 'time': 464 logger.warn("Specifying 'time' as a variable is dangerous! Are you " 465 "sure you know what you're doing?") 466 elif id == 'default': 467 logger.warn("'default' is a reserved keyword in C. This will cause " 468 "problems using the C-based integrator.") 469 elif id[0].isdigit(): 470 raise ValueError("The id %s is invalid. ids must not start with a " 471 "number." % id) 472 if id in self.variables.keys()\ 473 or id in self.reactions.keys()\ 474 or id in self.functionDefinitions.keys()\ 475 or id in self.events.keys()\ 476 or id in self.constraints.keys()\ 477 or id == self.id: 478 raise ValueError, ('The id %s is already in use!' % id)

479

480 - def set_id(self, id):

481 """ 482 Set the id of this Network 483 """ 484 if id != self.id: 485 self._checkIdUniqueness(id) 486 self.id = id

487

488 - def get_id(self):

489 """ 490 Get the id of this Network. 491 """ 492 return self.id

493

494 - def set_name(self, name):

495 """ 496 Set the name of this Network 497 """ 498 self.name = name

499

500 - def get_name(self):

501 """ 502 Get the name of this Network. 503 """ 504 return self.name

505

506 - def set_deterministic(self):

507 """ 508 Disables the stochastic simulation of network dynamics 509 """ 510 if hasattr(self, 'stochastic'): delattr(self, 'stochastic')

511

512 - def set_stochastic(self, seed=None, fill_dt=None, rmsd=None):

513 """ 514 set_stochastic enables the stochastic simulation of the network 515 dynamics (instead of deterministic integration). This will discretize 516 the dynamic variables in the model (making them integers). This 517 simulation has only limited support for events*, and does not 518 implement algebraic rules** or rates rules**. We use a kinetic MC 519 algorithm (aka Gillespie algorithm). The complementary function 520 'set_deterministic' may be used to disable stochastic simulations of 521 network dynamics. Please note the results returned for exactly the 522 times asked are interpolations between the surrounding points, thus 523 concentrations may be floats instead of integers. 524 525 * Since there is no exact way to implement events in a stochastic 526 simulation, only a very primative method is present. This algorithm 527 currently only supports time triggered events (more complex triggers 528 may be added in the future) and fires them as soon after that time 529 is passed as possible. Please note that it is the users 530 responsibility to ensure the events are firing in a sensible 531 way. Non-sensible output will arise when the time steps are too 532 large and events are not fired near their target trigger times. The 533 fired times are logged in the logger at the debug level. 534 535 ** Possible future implementation, although parsing rate rules have the 536 possible issue of creating negative reaction rates (which have no 537 physical meaning). Add reactions instead of rate rules for stochastic 538 simulations. 539 540 Inputs: 541 (seed=None, fill_dt=None, rmsd=None) 542 seed -- RNG seed used in the simulation. If none, generated from the 543 system time and os.urandom (if available). 544 fill_dt -- Fill in the trajectory at the requested time interval. This 545 is the simplest way to fill in a trajectory. 546 rmsd -- Maximum root mean squared distance in the dynamic variable 547 space before taking a data point (a better trajectory filler). 548 549 Outputs: None 550 """ 551 if seed==None: 552 def rehash(key): # From numpy's mtrand module 553 key += ~(key << 15) 554 key ^= (key >> 10) 555 key += (key << 3) 556 key ^= (key >> 6) 557 key += ~(key << 11) 558 key ^= (key >> 16) 559 return key

560 561 t = time.time() 562 ts = int(t) 563 tms = int((t-float(ts))*1000000.0) 564 seed = rehash(ts)^rehash(tms) 565 566 try: 567 import struct 568 seed ^= rehash( 569 struct.unpack('i', os.urandom(struct.calcsize('i')))[0]) 570 except: pass 571 572 seed = int(seed % sys.maxint) 573 574 if rmsd==None: 575 rmsd = _double_max_ 576 577 self.stochastic = {'seed':seed, 'reseed':True, 578 'fill_dt':fill_dt, 'rmsd':rmsd} 579 580 if hasattr(self, 'periodic'): 581 logger.warn('Cannot find periodic solutions with a stochastic solution, disabling periodic solution finding') 582 delattr(self, 'periodic') 583 self.remove_component('_tau')

584

585 - def set_periodic(self, period, xtol=0.001, maxfun=100, phase=None, 586 minVel=None, level=2, log=False, rel=False):

587 """ 588 set_periodic ensures that a period solution has satisfactorily 589 approached a limit cycle before returning the solution. Please note that 590 this has not been implemented for sensitivity calculations. 591 592 Inputs: 593 (period, xtol=0.001, maxfun=100, phase=None, minVel=None, level=2, 594 log=False, rel=False) 595 period -- initial guess of the period (required) 596 xtol -- tolerance required for convergence (all chemicals). If set too 597 high, numerical (in)precision may cause no stable limit cycle. 598 maxfun -- maximum number of iterations (one trajectory found in each) 599 phase -- set to a ('name', value) tuple to set the t=0 solution to a 600 particular state in systems which no forcing. Value may be 601 'max','min', or a float concentration (careful with the later). 602 minVel -- used to detect fixed points: 603 when the (max IC velocity)**2 < minVel the integration 604 is halted, the stability set to True, and the period set to 605 0.0 (indicating nonperiodic, fixed point soln). 606 level -- correlation level (>=1) to use to find the period. For example, 607 if level=2, the difference between x(0)-x(1*tau) and 608 x(0)-x(2*tau) is used. 609 log -- Set to True to use logarithms of state variables, constraining 610 the solution to non-zero values 611 rel -- Set to True to use relative differences (x1-x0)/min(x0,x1) in 612 the limit cycle period search. 613 614 Outputs: (stored in net.periodic dictionary) 615 tol -- tolerance of the solution (rmsd of ic1 to ic0 for each level) 616 period -- actual period of oscillation 617 stable -- whether a stable solution was found 618 feval -- number of iterations 619 stableIC -- stable initial condition 620 621 Additionally, the period is accessible during the calculation and 622 afterwards via the '_tau' variable (useful for assignments which 623 rely on the period). 624 """ 625 self.periodic = {'xtol': xtol, 'maxfun':maxfun, 'phase':phase, 626 'level':int(level), 'log':log, 'rel':rel, 627 'feval': 0, 'period': period, 'minVel': minVel, 628 'stable': False, 'stableIC': {}} 629 if hasattr(self, 'stochastic'): 630 logger.warn('Cannot find periodic solutions with a stochastic solution, switching to dynamic solution') 631 self.set_deterministic() 632 self.add_parameter('_tau', period, is_constant=True, 633 is_optimizable=False)

634

635 - def _iter_limit_cycle(self, params, varNames, s0):

636 """ 637 Internal function used to integrate a trajectory and find the point 638 nearest the initial point (defined as the root sum squared distance) 639 in the period range of [0.475*tau, 1.5*tau]. The point found and period 640 (tau) is returned. 641 642 Inputs: 643 (params, varNames, s0) 644 params -- Network parameters 645 varNames -- Names corresponding to the values in ln_x0 646 s0 -- Initial state (period and conditions) where s0[0]=period and 647 s0[1:]=initial conditions 648 649 Outputs: 650 (s1) 651 s1 -- New state found (period and conditions defined same as input) 652 """ 653 if self.periodic['rel']: 654 rmsd = lambda _s1, _s2: \ 655 scipy.sqrt(scipy.average(((_s1-_s2) / \ 656 scipy.minimum(_s1,_s2))**2)) 657 elif self.periodic['log']: 658 rmsd = lambda _s1, _s2: \ 659 scipy.sqrt(scipy.average((scipy.log(_s1/_s2)**2))) 660 else: 661 rmsd = lambda _s1, _s2: scipy.sqrt(scipy.average((_s1-_s2)**2)) 662 663 # Set the initial condition 664 s0 = scipy.array(s0) 665 tau0, x0 = s0[0], s0[1:] 666 for name, value in zip(varNames, x0): 667 self.set_var_ic(name, value) 668 669 # Set the period search interval, restrictions: 670 # 1) first guess by fminbound is at tau0 671 # 2) Lower bound < 0.5*tau0 (to possibly detect double periods) 672 # 3) Upper bound < 2.0*tau0 (to prevent double periods) 673 a = 0.475 674 b = a + (1.0-a)/(0.5*(3.0-scipy.sqrt(5.0))) # ~1.85 675 676 # Find the trajectory 677 max_time = b * tau0 * float(self.periodic['level']) 678 self.trajectory = self.integrate([0.,max_time], params, addTimes = True) 679 self.trajectory.build_interpolated_traj() 680 ev_inter_traj = self.trajectory.evaluate_interpolated_traj 681 682 # Find the period between 1/2*tau0 and 3/2*tau0 683 # Without the squeeze, this failed on some versions of numpy when 684 # t = array([0.2]) 685 s = lambda t: scipy.array([scipy.squeeze(t)] + \ 686 [ev_inter_traj(name, t) for name in varNames]) 687 688 def __iter_cost(t): 689 t = scipy.squeeze(abs(t)) 690 res = scipy.array([t], scipy.float_) 691 for lvl in range(1, self.periodic['level']+1): 692 res = scipy.concatenate((res, s(t*float(lvl))[1:])) 693 res0 = scipy.concatenate(([tau0], list(x0)*self.periodic['level'])) 694 return rmsd(res, res0)

695 696 tau1 = scipy.optimize.fminbound(__iter_cost, a * tau0, b * tau0) 697 tau1 = abs(tau1) 698 s1 = s(tau1) 699 700 # Update the periodic dictionary 701 self.periodic['feval'] += 1 702 self.periodic['period'] = tau1 703 self.set_var_val('_tau', tau1) 704 self.periodic['tol'] = __iter_cost(tau1) 705 if self.periodic['tol'] < self.periodic['xtol']: 706 self.periodic['stable']=True 707 708 logger.debug('Limit Cycle Iteration #%i: period = %f, tol = %f, %s'% \ 709 (self.periodic['feval'], self.periodic['period'], 710 self.periodic['tol'], 711 '%sstable'%('un'*(not self.periodic['stable'])))) 712 713 return s1 714

715 - def _eliminate_slowest_mode(self, s0, s1, s2):

716 """ 717 Internal function which finds and eliminates the slowest decaying 718 mode in the approach to a limit cycle: 719 si = [tau_i, x0_0, .., x0_i] 720 F0 = s_star + vector_1 + ... 721 F1 = s_star + lambda_1 * vector_1 + ... 722 F2 = s_star + lambda_1^2 * vector_1 + ... 723 724 lambda_1 = (F2-F1)*(F1-F0)/|F1-F0|^2 725 vector_1 = (F2-F1)/(lambda_1 * (lambda_1 - 1.0)) 726 727 Inputs: (s0, s1, s2) 728 Three vectors representing three consecutive steps towards a limit cycle 729 730 Ouptuts: (s_star) 731 s_star = s2 - lambda_1^2 * vector_1 732 """ 733 x0, x1, x2 = s0[1:], s1[1:], s2[1:] 734 lambda_1 = sum( (x2-x1) * (x1-x0) ) / sum( (x1-x0)**2 ) 735 vector_1 = (x2-x1)/(lambda_1 * (lambda_1-1.)) 736 s = scipy.array(s2) # make a copy of s2 737 s[1:] -= lambda_1*lambda_1*vector_1 738 return s

739

740 - def _find_limit_cycle(self, params):

741 """ 742 Internal function to find a stable limit cycle (if one exists), given a 743 parameter set, an estimated period, and a set of initial conditions. 744 This function calculates two iterations (which are hopefully approaching 745 a limit cycle), and use this to estimate a point on the limit cycle - an 746 estimation which is kept if it is closer to the limit cycle. 747 748 Inputs: 749 (params,) 750 params -- model parameters 751 752 Outputs: None 753 """ 754 self.periodic['stable'] = False 755 self.periodic['feval'] = 0 756 757 # We only want dynamic variables whose value is not assigned 758 varNames = [name for name in self.dynamicVars.keys() 759 if name in self.species.keys()] 760 761 # The initial period and state 762 tau0 = self.periodic['period'] 763 x0 = [self.periodic['stableIC'].get(name, self.get_var_ic(name)) 764 for name in varNames] 765 s0 = scipy.array([tau0]+x0) 766 767 while not self.periodic['stable'] and \ 768 self.periodic['feval'] < self.periodic['maxfun']: 769 770 # Perform the first two searched 771 s1 = self._iter_limit_cycle(params, varNames, s0) 772 if self.periodic['stable']: 773 s0 = s1 774 break 775 776 s2 = self._iter_limit_cycle(params, varNames, s1) 777 tol2 = self.periodic['tol'] 778 if self.periodic['stable']: 779 s0 = s2 780 break 781 782 if self.periodic['log']: 783 ls0, ls1, ls2 = scipy.log(s0), scipy.log(s1), scipy.log(s2) 784 s_star = scipy.exp(self._eliminate_slowest_mode(ls0, ls1, ls2)) 785 else: 786 s_star = self._eliminate_slowest_mode(s0, s1, s2) 787 788 # Predict a point on the orbit, and choose the best point to 789 # use as the next round (either the last search #2 or the 790 # prediction) 791 try: 792 s3 = self._iter_limit_cycle(params, varNames, s_star) 793 794 # Did we find a better solution? 795 if self.periodic['tol'] < tol2: s0 = s3 796 else: s0 = s2 # Prediction (s_star) was worse 797 except: 798 s0 = s2 799 800 # Found NANs, cut losses and run 801 if scipy.any(scipy.isnan(s0)): 802 self.periodic['stable'] = False 803 break 804 805 # Test for a fixed point, stopping if one is found 806 if self.periodic['minVel'] is not None: 807 for name, val in zip(varNames, s0[1:]): 808 self.set_var_ic(name, val) 809 810 maxVel = max( 811 [iV**2 for name, iV in zip(self.dynamicVars.keys(), 812 self.get_initial_velocities()) 813 if name in varNames]) 814 815 if maxVel < self.periodic['minVel']: 816 self.periodic['period']=0.0 817 self.set_var_val('_tau', 0.0) 818 self.periodic['stable']=True 819 break 820 821 822 if self.periodic['stable']: 823 # Record the initial conditions (if stable) 824 for name, value in zip(varNames, s0[1:]): 825 self.periodic['stableIC'][name] = value 826 self.set_var_ic(name, value) 827 828 # Set the phase of the cycle (or else return, we're done) 829 if self.periodic['phase'] is not None: 830 # Define the cost 831 pN, pV = self.periodic['phase'] 832 _eit = self.trajectory.evaluate_interpolated_traj 833 834 if pV=='min': cost = lambda t: _eit(pN, t) 835 elif pV=='max': cost = lambda t: 1.0 / _eit(pN, t) 836 else: cost = lambda t: (pV - _eit(pN, t))**2 837 838 # Find the best-fitting point 839 tau = self.periodic['period'] 840 lowerB = float(self.periodic['level'] - 1) * tau 841 upperB = float(self.periodic['level']) * tau 842 t0 = scipy.optimize.fminbound(cost, lowerB, upperB) 843 844 # Set the initial point, fall through and integrate 845 for name in varNames: 846 self.periodic['stableIC'][name] = _eit(name, t0) 847 self.set_var_ic(name, _eit(name, t0))

848 849 850 # 851 # Methods to become a 'SloppyCell.Model' 852 #

853 - def calculate(self, vars, params=None):

854 self.Calculate(vars, params) 855 return self.GetResult(vars)

856

857 - def Calculate(self, vars, params = None):

858 # Add in the times required by all the variables, then convert back to 859 # a sorted list. 860 # Make sure we start from t = 0 861 t = sets.Set([0]) 862 ret_full_traj=False 863 for var, times in vars.items(): 864 if var == 'full trajectory': 865 ret_full_traj = True 866 t.union_update(sets.Set(times)) 867 elif var.endswith('_maximum') or var.endswith('_minimum'): 868 t1,t2 = times 869 if t1 is not None: 870 t.add(t1) 871 if t2 is not None: 872 t.add(t2) 873 elif self.variables.has_key(var): 874 t.union_update(sets.Set(times)) 875 else: 876 raise ValueError('Unknown variable %s requested from network %s' 877 % (var, self.get_id())) 878 879 # This takes care of the normal data points 880 t = list(t) 881 t.sort() 882 883 if hasattr(self, 'periodic'): 884 self._find_limit_cycle(params) 885 if self.periodic['period'] > max(t): 886 t.append(self.periodic['period']) 887 elif hasattr(self, 'stochastic'): 888 self.trajectory = self.integrateStochastic(t, params=params) 889 return 890 891 self.trajectory = self.integrate(t, params=params, 892 addTimes=ret_full_traj)

893

894 - def CalculateSensitivity(self, vars, params):

895 t = sets.Set([0]) 896 897 for var,times in vars.items(): 898 t.union_update(sets.Set(times)) 899 900 t = list(t) 901 t.sort() 902 903 self.ddv_dpTrajectory = self.integrateSensitivity(t,params=params, 904 addTimes=True) 905 self.trajectory = self.ddv_dpTrajectory

906

907 - def GetName(self):

908 return self.get_id()

909

910 - def GetParameters(self):

911 return KeyedList([(var.id, var.value) for var in 912 self.optimizableVars.values()])

913

914 - def GetParameterTypicalValues(self):

915 return KeyedList([(var.id, var.typicalValue) for var in 916 self.optimizableVars.values()])

917

918 - def GetResult(self, vars):

919 result = {} 920 times = self.trajectory.timepoints 921 for id in vars: 922 if id == 'full trajectory': 923 self.trajectory.build_interpolated_traj() 924 result[id] = self.trajectory 925 if id.endswith('_maximum') or id.endswith('_minimum'): 926 # This is an extremum search. 927 # Which variable are we looking for? 928 var = '_'.join(id.split('_')[:-1]) 929 minTime, maxTime = vars[id] 930 931 # First, check though the trajectory, to see if the extremum 932 # is in there. 933 t = self.trajectory.get_times() 934 # Restrict timespan in which to search. 935 if minTime is not None: 936 min_t_ii = t.searchsorted(minTime) 937 else: 938 min_t_ii = 0 939 if maxTime is not None: 940 max_t_ii = t.searchsorted(maxTime)-1 941 else: 942 max_t_ii = len(t) 943 vartraj = self.trajectory.get_var_traj(var) 944 if id.endswith('_maximum'): 945 var_val_ii = vartraj[min_t_ii:max_t_ii+1].argmax()\ 946 + min_t_ii 947 elif id.endswith('_minimum'): 948 var_val_ii = vartraj[min_t_ii:max_t_ii+1].argmin()\ 949 + min_t_ii 950 var_val = vartraj[var_val_ii] 951 t_val = t[var_val_ii] 952 953 # Now check through the events. If we're running under 954 # full_speed the trajectory can be very sparse. If the use has 955 # set up events to tag potential extrema, this will use those. 956 curr_vals = self.get_var_vals() 957 event_info = self.trajectory.event_info 958 for ii in range(len(event_info[0])): 959 time,dynVarVals = event_info[0][ii], event_info[1][ii] 960 self.updateVariablesFromDynamicVars(dynVarVals, time) 961 eval = self.get_var_val(var) 962 if id.endswith('_maximum') and eval > var_val\ 963 and time >= minTime and time <= maxTime: 964 var_val = eval 965 t_val = time 966 elif id.endswith('_minimum') and eval < var_val\ 967 and time >= minTime and time <= maxTime: 968 var_val = eval 969 t_val = time 970 self.set_var_vals(curr_vals) 971 972 result[id] = {(minTime, maxTime):(t_val,var_val)} 973 elif self.variables.has_key(id): 974 traj = self.trajectory.getVariableTrajectory(id) 975 result[id] = dict(zip(times, traj)) 976 977 return result

978

979 - def GetSensitivityResult(self, vars):

980 # note: returns all the timepoints we have, not just 981 # the requested ones (which should be a subset of all 982 # the timepoints) 983 result = {} 984 times = self.ddv_dpTrajectory.get_times() 985 for id in vars.keys(): 986 result[id] = {} 987 for tIndex, t in enumerate(times): 988 result[id][t] = {} 989 for optparam in self.optimizableVars.keys(): 990 result[id][t][optparam] = \ 991 self.ddv_dpTrajectory.get_var_traj((id,optparam))[tIndex] 992 return result

993

994 - def integrate(self, times, params=None, returnEvents=False, addTimes=True, 995 rtol = None):

996 if self.add_tail_times or addTimes: 997 times = scipy.concatenate((times, [1.05*times[-1]])) 998 return Dynamics.integrate(self, times, params=params, 999 fill_traj=(self.add_int_times or addTimes))

1000

1001 - def integrateSensitivity(self, times, params = None, 1002 returnEvents = False, addTimes = True, 1003 rtol=None):

1004 if self.add_tail_times: 1005 times = scipy.concatenate((times, [1.05*times[-1]])) 1006 return Dynamics.integrate_sensitivity(self, times, params, rtol, 1007 fill_traj=self.add_int_times)

1008

1009 - def integrateStochastic(self, times, params=None):

1010 if self.stochastic['fill_dt'] is not None: 1011 times = sets.Set(times) 1012 times.union_update(scipy.arange(min(times), max(times), 1013 self.stochastic['fill_dt'])) 1014 times = list(times) 1015 times.sort() 1016 1017 if times[0]==0.: # Only reset if specifically requested (not by event) 1018 self.resetDynamicVariables() 1019 1020 # Add in the event times (only if they don't extend the trajectoy!) 1021 times = sets.Set(times) 1022 times.union_update([_.triggeringTime for _ in self.events 1023 if _.triggeringTime<max(times)]) 1024 times = list(times) 1025 times.sort() 1026 1027 if times[0]==0.: 1028 times = times[1:] 1029 1030 if params is not None: self.update_optimizable_vars(params) 1031 1032 self.compile() 1033 1034 # Test that we have a working integrator 1035 body = self._dynamic_funcs_python.get('integrate_stochastic_tidbit') 1036 if 'pass' in body.split('\n')[-1]: 1037 err_body = '\n'.join( 1038 ['Integrate stochastic failed due to the problem(s):']+ 1039 body.split('\n')[1:-1]) 1040 raise RuntimeError(err_body) 1041 1042 dv=scipy.array([self.get_var_val(_) for _ in self.dynamicVars.keys()]) 1043 cv = self.constantVarValues 1044 1045 trajectory = Trajectory_mod.Trajectory(self, holds_dt=0, const_vals=cv) 1046 trajectory.appendFromODEINT(scipy.array([0.0]), scipy.array([dv])) 1047 1048 # VERY simple event handling... 1049 events_occurred = [] 1050 pendingEvents = {} 1051 for event in self.events.values(): 1052 pendingEvents.setdefault(event.triggeringTime, []) 1053 pendingEvents[event.triggeringTime].append(event) 1054 1055 t, tInd = 0.0, 0 1056 while t < times[-1]: 1057 # Get the next chunk 1058 t, dv, tout, dvout = self.integrate_stochastic_tidbit( 1059 self.stochastic['seed'], self.stochastic['reseed'], 1060 t, dv, cv, self.stochastic['rmsd'], times[tInd]) 1061 self.stochastic['reseed']=False 1062 1063 # Fill in trajectory and execute any events that occurred 1064 while tInd<len(times) and t >= times[tInd]: 1065 # Extrapolate the dynamic vars, if not the stop_time requested 1066 if times[tInd]!=tout: 1067 dt = t - trajectory.timepoints[-1] 1068 dvout = trajectory.values[-1]+(dv-trajectory.values[-1])/\ 1069 dt*(times[tInd] - trajectory.timepoints[-1]) 1070 1071 trajectory.appendFromODEINT(scipy.array([times[tInd]]), 1072 scipy.array([dvout])) 1073 1074 # Is this where an event fires? 1075 if times[tInd] in pendingEvents.keys(): 1076 events = pendingEvents.pop(times[tInd]) 1077 for event in events: 1078 # Execute twice b/c we want to update the dvout at 1079 # the fired time in addition to updating the dv for 1080 # any future extrapolations. This is really a product 1081 # of not being able to ensure we haven't overstepped 1082 # our desired timepoint and crossed the next. 1083 dvout = self.executeEvent(event, times[tInd], dvout, \ 1084 dvout, times[tInd]) 1085 dv = self.executeEvent(event, times[tInd], dvout, dv, t) 1086 events_occurred.append(event) 1087 logger.debug('Executed event %s in net %s at t=%s'%\ 1088 (event.id, self.id, t)) 1089 trajectory.appendFromODEINT(scipy.array([times[tInd]]), 1090 scipy.array([dvout])) 1091 1092 tInd += 1 1093 1094 if trajectory.timepoints[-1]<t: 1095 trajectory.appendFromODEINT(scipy.array([t]), 1096 scipy.array([dv])) 1097 1098 return trajectory

1099

1100 - def _sub_for_piecewise(self, expr, time):

1101 """ 1102 Runs through expr and replaces all piecewise expressions by the 1103 clause that is currently active. 1104 """ 1105 if not expr.count('piecewise('): 1106 return expr 1107 funcs_used = ExprManip.extract_funcs(expr) 1108 for func, args in funcs_used: 1109 if func == 'piecewise': 1110 ast = ExprManip.strip_parse(expr) 1111 ast = self._sub_for_piecewise_ast(ast, time) 1112 return ExprManip.ast2str(ast) 1113 return expr

1114

1115 - def _sub_for_piecewise_ast(self, ast, time):

1116 if isinstance(ast, ExprManip.AST.CallFunc)\ 1117 and ExprManip.ast2str(ast.node) == 'piecewise': 1118 # If our ast is a piecewise function 1119 conditions = [cond for cond in ast.args[1::2]] 1120 clauses = [clause for clause in ast.args[:-1:2]] 1121 if len(ast.args)%2 == 1: 1122 # odd # of arguments implies otherwise clause 1123 otherwise = ast.args[-1] 1124 else: 1125 otherwise = None 1126 1127 for cond_ast, expr_ast in zip(conditions, clauses): 1128 # For each condition, check whether or not it is True. 1129 # If so, return the corresponding clause. 1130 subbed_cond = self._sub_for_piecewise_ast(cond_ast, time=time) 1131 cond_str = ExprManip.ast2str(subbed_cond) 1132 if self.evaluate_expr(cond_str, time): 1133 return self._sub_for_piecewise_ast(expr_ast, time=time) 1134 else: 1135 # If none of our conditions were True, return the otherwise 1136 # clause if we have one. 1137 if otherwise is not None: 1138 return self._sub_for_piecewise_ast(otherwise, time=time) 1139 else: 1140 raise ValueError('No True condition and no otherwise ' 1141 " clause in '%s'." 1142 % ExprManip.ast2str(ast)) 1143 return ExprManip.AST.recurse_down_tree(ast, self._sub_for_piecewise_ast, 1144 (time,))

1145

1146 - def evaluate_expr(self, expr, time=0, var_vals=None):

1147 """ 1148 Evaluate the given expression using the current values of the network 1149 variables. 1150 """ 1151 try: 1152 return float(expr) 1153 except ValueError: 1154 # We create a local_namespace to evaluate the expression in that 1155 # maps variable ids to their current values 1156 expr = self._sub_for_piecewise(expr, time) 1157 if var_vals is None: 1158 vars_used = ExprManip.extract_vars(expr) 1159 var_vals = [(id, self.get_var_val(id)) for id in vars_used 1160 if id != 'time'] 1161 var_vals = dict(var_vals) 1162 var_vals['time'] = time 1163 local_namespace = var_vals 1164 local_namespace.update(self.namespace) 1165 # We strip whitespace, just for convenience 1166 return eval(expr.strip(), local_namespace, {})

1167 1168 # 1169 # Methods to get and set object properties 1170 #

1171 - def set_var_typical_val(self, id, value):

1172 """ 1173 Set the typical value for a variable. 1174 """ 1175 var = self.get_variable(id) 1176 var.typicalValue = value

1177

1178 - def get_var_typical_val(self, id):

1179 """ 1180 Return the typical value for a variable. 1181 """ 1182 return self.get_variable(id).typicalValue

1183

1184 - def get_var_typical_vals(self):

1185 """ 1186 Return the variable typical values as a KeyedList. 1187 """ 1188 return KeyedList([(id, self.get_var_typical_val(id)) for id in 1189 self.variables.keys()])

1190

1191 - def set_var_ic(self, id, value, warn=True, update_constants=True):

1192 """ 1193 Set the initial condition of the variable with the given id. 1194 """ 1195 if warn and id in self.assignedVars.keys(): 1196 logger.warn('WARNING! Attempt to assign an initial condition to ' 1197 'the variable %s, which is determined by an assignment ' 1198 'rule. This is a meaningless operation. Instead, ' 1199 'change the initial condition of one or more of the ' 1200 'components in the rule: %s' % 1201 (id, self.assignmentRules.get(id))) 1202 1203 var = self.get_variable(id) 1204 var.initialValue = value 1205 if var.is_constant: 1206 var.value = value 1207 if update_constants: 1208 self.constantVarValues = [self.evaluate_expr(var.value) for var 1209 in self.constantVars.values()] 1210 self.constantVarValues = scipy.array(self.constantVarValues)

1211

1212 - def get_var_ic(self, id):

1213 """ 1214 Return the initial condition for a variable 1215 """ 1216 return self.get_variable(id).initialValue

1217

1218 - def get_var_ics(self):

1219 """ 1220 Return the variable initial conditions as a KeyedList 1221 """ 1222 return KeyedList([(id, self.get_var_ic(id)) for id in 1223 self.variables.keys()])

1224

1225 - def set_var_vals(self, kl, time = 0):

1226 """ 1227 Set current variable values from a KeyedList or dictionary. 1228 """ 1229 for id, value in kl.items(): 1230 self.set_var_val(id, value, time, warn=False, do_assignments=False) 1231 self.updateAssignedVars(time)

1232

1233 - def get_var_val(self, id):

1234 """ 1235 Return the current value of a variable 1236 """ 1237 val = self.get_variable(id).value 1238 return self.evaluate_expr(val)

1239

1240 - def get_var_vals(self, ids=None):

1241 """ 1242 Return the current variable values as a KeyedList 1243 1244 ids -- List o variable ids to return values for. If None, all variables 1245 are used. 1246 """ 1247 if ids is None: 1248 ids = self.variables.keys() 1249 return KeyedList([(id, self.get_var_val(id)) for id in ids])

1250

1251 - def get_initial_velocities(self):

1252 """ 1253 Returns the vector field evaluated at the initial conditions 1254 """ 1255 ics = [self.evaluate_expr(self.getInitialVariableValue(dvid)) 1256 for dvid in self.dynamicVars.keys()] 1257 # the evaluate_expr is in case an initial condition is a parameter 1258 y0, initialv = Dynamics.find_ics(y=ics, yp=scipy.ones(len(ics)), 1259 time=0, 1260 var_types=self._dynamic_var_algebraic, 1261 rtol=[1e-6]*len(ics), 1262 atol=[1e-6]*len(ics), 1263 constants=self.constantVarValues, 1264 net=self) 1265 return initialv

1266

1267 - def set_var_ics(self, kl):

1268 """ 1269 Set variable initial conditions from a KeyedList or dictionary. 1270 """ 1271 for id, value in kl.items(): 1272 if id in self.variables.keys(): 1273 self.set_var_ic(id, value, warn=False)

1274

1275 - def set_var_val(self, id, val, time=0, warn=True, do_assignments=True):

1276 """ 1277 Set the current stored value of the variable with the given id. 1278 """ 1279 1280 if warn and self.assignedVars.has_key(id): 1281 logger.warn('WARNING! Attempt to assign a value to the variable ' 1282 '%s, which is determined by an assignment rule. This ' 1283 'is a meaningless operation. Instead, change the value ' 1284 'of one or more of the components in the rule: %s' 1285 % (id, self.assignmentRules.get(id))) 1286 1287 var = self.get_variable(id) 1288 var.value = val 1289 if do_assignments: 1290 self.updateAssignedVars(time)

1291 1292 setTypicalVariableValue = set_var_typical_val 1293

1294 - def get_variable(self, id):

1295 """ 1296 Return the class instance with a given variable id. 1297 """ 1298 var = self.variables.get(id) 1299 if var: 1300 return var 1301 else: 1302 raise KeyError('Variable %s not found in network %s!' 1303 % (id, self.get_id()))

1304

1305 - def update_optimizable_vars(self, params):

1306 """ 1307 Update the net's optimizable vars from a passed-in KeyedList, 1308 dictionary, or sequence. 1309 1310 Only those variables that intersect between the net and params are 1311 changed if a KeyedList or dict is passed in. 1312 """ 1313 if hasattr(params, 'get'): 1314 inBoth = sets.Set(self.optimizableVars.keys()) 1315 inBoth = inBoth.intersection(sets.Set(params.keys())) 1316 for id in inBoth: 1317 self.set_var_ic(id, params.get(id), update_constants=False) 1318 elif len(params) == len(self.optimizableVars): 1319 for ii, id in enumerate(self.optimizableVars.keys()): 1320 self.set_var_ic(id, params[ii], update_constants=False) 1321 else: 1322 raise ValueError('Passed in parameter set does not have the proper ' 1323 'length!') 1324 1325 self.constantVarValues = [self.evaluate_expr(var.value) for var in 1326 self.constantVars.values()] 1327 self.constantVarValues = scipy.array(self.constantVarValues)

1328 1329 getInitialVariableValue = get_var_ic 1330

1331 - def getDynamicVarValues(self):

1332 # We need to evaluate_expr here to handle ones that are assigned to 1333 # parameters 1334 return scipy.array([self.evaluate_expr(var.value) 1335 for var in self.dynamicVars.values()])

1336

1337 - def resetDynamicVariables(self):

1338 # Resets all dynamical variables to their initial values. This is 1339 # a little complex because the initial value may be a function of 1340 # other values. Thus we skip those ones on the first pass through. 1341 pending = [] 1342 for id, var in self.dynamicVars.items(): 1343 if isinstance(var.initialValue, str): 1344 pending.append((id, var)) 1345 else: 1346 var.value = var.initialValue 1347 1348 # We need to update the assigned variables both before and after 1349 # we do the pending dynamic vars. 1350 self.updateAssignedVars(time = 0) 1351 1352 for id, var in pending: 1353 self.dynamicVars.getByKey(id).value = \ 1354 self.evaluate_expr(var.initialValue, 0) 1355 1356 self.updateAssignedVars(time = 0)

1357

1358 - def set_dyn_var_ics(self, values):

1359 for ii, id in enumerate(self.dynamicVars.keys()): 1360 if hasattr(values, 'get'): 1361 self.set_var_ic(id, values.get(id)) 1362 else: 1363 self.set_var_ic(id, values[ii])

1364

1365 - def updateVariablesFromDynamicVars(self, values, time):

1366 for ii in range(len(self.dynamicVars)): 1367 self.dynamicVars[ii].value = values[ii] 1368 self.updateAssignedVars(time)

1369

1370 - def updateAssignedVars(self, time):

1371 to_get = sets.Set(self.variables.keys()) 1372 to_get.difference_update(sets.Set(self.assignedVars.keys())) 1373 var_vals = [(id, self.get_var_val(id)) for id in to_get] 1374 var_vals = dict(var_vals) 1375 var_vals['time'] = time 1376 for id, rhs in self.assignmentRules.items(): 1377 assigned_val = self.evaluate_expr(rhs, time, var_vals=var_vals) 1378 self.assignedVars.getByKey(id).value = assigned_val 1379 var_vals[id] = assigned_val

1380

1381 - def set_var_optimizable(self, id, is_optimizable):

1382 self.get_variable(id).is_optimizable = is_optimizable 1383 self._makeCrossReferences()

1384

1385 - def set_var_constant(self, id, is_constant):

1386 self.get_variable(id).is_constant = is_constant 1387 self._makeCrossReferences()

1388

1389 - def get_var_constant(self, id):

1390 return self.get_variable(id).is_constant

1391 1392 # 1393 # Generate the differential equations and functions to calculate them. 1394 # 1395

1396 - def _makeDiffEqRHS(self):

1397 logger.debug('Making diff equation rhs') 1398 diff_eq_terms = {} 1399 1400 for rxn_id, rxn in self.reactions.items(): 1401 rateExpr = rxn.kineticLaw 1402 logger.debug('Parsing reaction %s.' % rxn_id) 1403 1404 for reactantId, dReactant in rxn.stoichiometry.items(): 1405 if self.get_variable(reactantId).is_boundary_condition or\ 1406 self.get_variable(reactantId).is_constant or\ 1407 self.assignmentRules.has_key(reactantId): 1408 # Variables that are boundary conditions, are constant, or 1409 # are assigned aren't modified by reactions, so we move on 1410 # to the next. 1411 continue 1412 1413 term = '(%s) * (%s)' % (dReactant, rateExpr) 1414 diff_eq_terms.setdefault(reactantId, []) 1415 diff_eq_terms[reactantId].append(term) 1416 1417 self.diff_eq_rhs = KeyedList() 1418 for id in self.dynamicVars.keys(): 1419 if self.rateRules.has_key(id): 1420 self.diff_eq_rhs.set(id, self.rateRules.get(id)) 1421 else: 1422 # We use .get to return a default of ['0'] 1423 rhs = '+'.join(diff_eq_terms.get(id, ['0'])) 1424 self.diff_eq_rhs.set(id, ExprManip.simplify_expr(rhs))

1425

1426 - def _make_res_function(self):

1427 py_body = [] 1428 py_body.append('def res_function(time, dynamicVars, yprime, ' 1429 'constants):') 1430 py_body.append('dynamicVars = scipy.asarray(dynamicVars)') 1431 py_body.append('yprime = scipy.asarray(yprime)') 1432 py_body.append('constants = scipy.asarray(constants)') 1433 py_body.append('') 1434 py_body.append('residual = scipy.empty(%i, scipy.float_)' 1435 % len(self.dynamicVars)) 1436 py_body.append('') 1437 self._add_assignments_to_function_body(py_body) 1438 py_body.append('') 1439 1440 c_body = [] 1441 # The C version needs to take a number of arguments that the f2py 1442 # wrapper will hide. 1443 # To be callable from Fortran it also needs 1) its name followed by an 1444 # underscore, and 2) to take all arguments as pointers. 1445 c_args = 'double *time_ptr, double *dynamicVars, double *yprime, '\ 1446 'double *cj_ptr, double *residual, int *ires_ptr, '\ 1447 'double *constants, int *ipar' 1448 c_body.append('void res_function_(%s){' % c_args) 1449 self._prototypes_c.set('res_function', 1450 'void res_function_(%s);' % c_args) 1451 c_body.append('double time = *time_ptr;') 1452 c_body.append('') 1453 self._add_assignments_to_function_body(c_body, in_c=True) 1454 c_body.append('') 1455 1456 # Make a list of algebraic rules for accessing later 1457 algebraicRuleList = self.algebraicRules.values() 1458 1459 # We keep a list of terms in our residual function for use building the 1460 # various derivative functions. 1461 self._residual_terms = [] 1462 # This list tells us whether a give dynamic variable is algebraic or not 1463 # It will be used in the integration. 1464 self._dynamic_var_algebraic = [] 1465 1466 # Loop over everything in the dynamicVars list 1467 for ii, (id, var) in enumerate(self.dynamicVars.items()): 1468 if self.algebraicVars.has_key(id): 1469 # It's an algebraic equation. Pop the first algebraic rule off 1470 # the list. 1471 rhs = algebraicRuleList.pop(0) 1472 py_body.append('# Residual function corresponding to an ' 1473 'algebraic equation') 1474 py_body.append('residual[%i] = %s' % (ii, rhs)) 1475 c_rhs = ExprManip.make_c_compatible(rhs) 1476 c_body.append('residual[%i] = %s;' % (ii, c_rhs)) 1477 self._residual_terms.append((rhs, None)) 1478 self._dynamic_var_algebraic.append(-1) 1479 else: 1480 rhs = self.diff_eq_rhs.getByKey(id) 1481 c_rhs = ExprManip.make_c_compatible(rhs) 1482 self._dynamic_var_algebraic.append(+1) 1483 py_body.append('# Residual function for %s' % id) 1484 if rhs != '0': 1485 py_body.append('residual[%i] = %s - yprime.item(%i)' % 1486 (ii, rhs, ii)) 1487 c_body.append('residual[%i] = %s - yprime[%i];' % 1488 (ii, c_rhs, ii)) 1489 self._residual_terms.append((rhs, id)) 1490 else: 1491 py_body.append('residual[%i] = -yprime.item(%i)' % (ii, ii)) 1492 c_body.append('residual[%i] = -yprime[%i];' % (ii, ii)) 1493 self._residual_terms.append((None, id)) 1494 1495 py_body.append('') 1496 py_body.append('return residual') 1497 py_body = '\n '.join(py_body) 1498 1499 c_body.append('}') 1500 c_body = os.linesep.join(c_body) 1501 1502 self._dynamic_funcs_python.set('res_function', py_body) 1503 self._dynamic_funcs_c.set('res_function', c_body) 1504 1505 return py_body

1506

1507 - def _make_root_func(self):

1508 len_root_func = 0 1509 1510 py_body = [] 1511 py_body.append('def root_func(time, dynamicVars, yprime, constants):') 1512 py_body.append('dynamicVars = scipy.asarray(dynamicVars)') 1513 py_body.append('yprime = scipy.asarray(yprime)') 1514 py_body.append('constants = scipy.asarray(constants)') 1515 py_body.append('') 1516 # We don't know the length of root_devs yet, so for now we'll 1517 # just insert a placeholder line. 1518 py_body.append('root_devs = scipy.empty(NEED_TO_FIX, scipy.float_)') 1519 py_body.append('') 1520 self._add_assignments_to_function_body(py_body, include_dts=True) 1521 py_body.append('') 1522 1523 c_body = [] 1524 c_args = 'int *neq_ptr, double *time_ptr, double *dynamicVars, '\ 1525 'double *yprime, int *nrt_ptr, double *root_devs, '\ 1526 'double *constants, int *ipar' 1527 c_body.append('void root_func_(%s){' % c_args) 1528 c_body.append('double time = *time_ptr;') 1529 c_body.append('') 1530 self._add_assignments_to_function_body(c_body, in_c=True, 1531 include_dts=True) 1532 self._prototypes_c.set('root_func', 1533 'void root_func_(%s);' % c_args) 1534 c_body.append('') 1535 1536 self.event_clauses = [] 1537 for ii, event in enumerate(self.events.values()+self.constraints.values()): 1538 trigger = event.trigger 1539 for func_name, func_vars, func_expr in self._logical_comp_func_defs: 1540 trigger = ExprManip.sub_for_func(trigger, func_name, func_vars, 1541 func_expr) 1542 1543 py_body.append('root_devs[%i] = (%s) - 0.5' % (len_root_func, 1544 trigger)) 1545 c_trigger = ExprManip.make_c_compatible(trigger) 1546 c_body.append('root_devs[%i] = (%s) - 0.5;' % (len_root_func, 1547 c_trigger)) 1548 self.event_clauses.append(trigger) 1549 len_root_func += 1 1550 1551 for ii, event in enumerate(self.events.values()+self.constraints.values()): 1552 event.sub_clause_indices = [] 1553 trigger = event.trigger 1554 for func_name, func_vars, func_expr in self._logical_comp_func_defs: 1555 trigger = ExprManip.sub_for_func(trigger, func_name, func_vars, 1556 func_expr) 1557 comparisons = ExprManip.extract_comps(trigger) 1558 if len(comparisons) > 1: 1559 for comp in comparisons: 1560 py_body.append('root_devs[%i] = (%s) - 0.5\n\t' 1561 % (len_root_func, comp)) 1562 c_comp = ExprManip.make_c_compatible(comp) 1563 c_body.append('root_devs[%i] = (%s) - 0.5;\n\t' 1564 % (len_root_func, c_comp)) 1565 self.event_clauses.append(comp) 1566 event.sub_clause_indices.append(len_root_func) 1567 len_root_func += 1 1568 1569 self.len_root_func = len_root_func 1570 # Now that we know how many roots we're looking for, go back and 1571 # insert the proper size for root_devs 1572 py_body[5] = ('root_devs = scipy.empty(%i, scipy.float_)' 1573 % len_root_func) 1574 py_body.append('') 1575 py_body.append('return root_devs') 1576 py_body = '\n\t'.join(py_body) 1577 1578 c_body.append('}') 1579 c_body = os.linesep.join(c_body) 1580 1581 self._dynamic_funcs_python.set('root_func', py_body) 1582 self._dynamic_funcs_c.set('root_func', c_body) 1583 1584 return py_body

1585

1586 - def _make_alg_deriv_func(self):

1587 # This function is used when we want to calculate consistent derivatives 1588 # for algebraic variables. It returns the derivative wrt time of 1589 # all the algebraic rules. Since the rules are <0=rule>, these should 1590 # all be zero when we have consistent values for the variable derivs. 1591 # E.g. If we have an algebraic rule 0 = f(x,y), then taking the time 1592 # derivative yields: 0 = df/dx*dx/dt + df/dy*dy/dt. This function 1593 # returns df/dx*dx/dt + df/dy*dy/dt and fsolve twiddles with the 1594 # dx/dt's that correspond to algebraic variables. It's a linear systems 1595 # of equations, so fsolve is kind of overkill, but it hasn't been a 1596 # slow point so far. 1597 py_body = [] 1598 py_body.append('def alg_deriv_func(alg_yp, dynamicVars, yp, time, ' 1599 'constants):') 1600 py_body.append('alg_yp = scipy.asarray(alg_yp)') 1601 py_body.append('dynamicVars = scipy.asarray(dynamicVars)') 1602 py_body.append('yp = scipy.asarray(yp)') 1603 py_body.append('constants = scipy.asarray(constants)') 1604 py_body.append('') 1605 py_body.append('alg_derivs_res = scipy.empty(%i, scipy.float_)' 1606 % len(self.algebraicRules)) 1607 py_body.append('') 1608 self._add_assignments_to_function_body(py_body) 1609 py_body.append('') 1610 1611 c_body = [] 1612 c_args = 'double *alg_yp, double *dynamicVars, double *yp, '\ 1613 'double *time_ptr, double *constants, double *alg_derivs_res' 1614 c_body.append('void alg_deriv_func_(%s){' % c_args) 1615 c_body.append('double time = *time_ptr;') 1616 c_body.append('') 1617 self._add_assignments_to_function_body(c_body, in_c=True) 1618 self._prototypes_c.set('alg_deriv_func', 1619 'void alg_deriv_func_(%s);' % c_args) 1620 c_body.append('') 1621 1622 for rule_ii, rule in enumerate(self.algebraicRules.values()): 1623 # Each rhs is a sum of drule/dA * dA/dt + drule/dB * dB/dt + ... 1624 rhs_terms = [] 1625 rhs_terms_c = [] 1626 rule_vars = ExprManip.extract_vars(rule) 1627 deriv_wrt_time = self.takeDerivative(rule, 'time', rule_vars) 1628 rhs_terms.append(deriv_wrt_time) 1629 rhs_terms_c.append(deriv_wrt_time) 1630 for dyn_id in self.dynamicVars.keys(): 1631 deriv = self.takeDerivative(rule, dyn_id, rule_vars) 1632 if deriv != '0': 1633 if self.algebraicVars.has_key(dyn_id): 1634 index = self.algebraicVars.index_by_key(dyn_id) 1635 rhs_terms.append('(%s)*alg_yp.item(%i)' % (deriv,index)) 1636 rhs_terms_c.append('(%s)*alg_yp[%i]' % (deriv,index)) 1637 else: 1638 index = self.dynamicVars.index_by_key(dyn_id) 1639 rhs_terms.append('(%s)*yp.item(%i)' % (deriv, index)) 1640 rhs_terms_c.append('(%s)*yp[%i]' % (deriv, index)) 1641 rhs = ' + '.join(rhs_terms) 1642 rhs_c = ' + '.join(rhs_terms_c) 1643 rhs_c = ExprManip.make_c_compatible(rhs_c) 1644 py_body.append('alg_derivs_res[%i] = %s' % (rule_ii, rhs)) 1645 c_body.append('alg_derivs_res[%i] = %s;' % (rule_ii, rhs_c)) 1646 1647 py_body.append('') 1648 py_body.append('return alg_derivs_res') 1649 py_body = '\n\t'.join(py_body) 1650 1651 c_body.append('}') 1652 c_body = os.linesep.join(c_body) 1653 1654 self._dynamic_funcs_python.set('alg_deriv_func', py_body) 1655 self._dynamic_funcs_c.set('alg_deriv_func', c_body) 1656 1657 return py_body

1658

1659 - def _make_alg_res_func(self):

1660 py_body = [] 1661 py_body.append('def alg_res_func(alg_vals, dynamicVars, ' 1662 'time, constants):') 1663 py_body.append('alg_vals = scipy.asarray(alg_vals)') 1664 py_body.append('dynamicVars = scipy.asarray(dynamicVars)') 1665 py_body.append('constants = scipy.asarray(constants)') 1666 py_body.append('') 1667 py_body.append('residual = scipy.zeros(%i, scipy.float_)' 1668 % len(self.algebraicVars)) 1669 1670 c_body = [] 1671 c_args = 'double *alg_vals, double *dynamicVars, double *time_ptr, '\ 1672 'double *constants, double *residual' 1673 self._prototypes_c.set('alg_res_func', 1674 'void alg_res_func_(%s);' % c_args) 1675 c_body.append('void alg_res_func_(%s){' % c_args) 1676 c_body.append('double time = *time_ptr;') 1677 c_body.append('') 1678 1679 # Copy our algebraic variable guesses into the appropriate slots of 1680 # dynamicVars 1681 for ii, key in enumerate(self.algebraicVars.keys()): 1682 var_index = self.dynamicVars.index_by_key(key) 1683 py_body.append('dynamicVars[%i] = alg_vals[%i]' 1684 % (var_index, ii)) 1685 c_body.append('dynamicVars[%i] = alg_vals[%i];' 1686 % (var_index, ii)) 1687 1688 # Now make all the assignments from dynamicVars to the variable ids 1689 py_body.append('') 1690 self._add_assignments_to_function_body(py_body) 1691 py_body.append('') 1692 c_body.append('') 1693 self._add_assignments_to_function_body(c_body, in_c=True) 1694 c_body.append('') 1695 1696 # Now evaluate all the algebraic rules. 1697 for ii, rhs in enumerate(self.algebraicRules.values()): 1698 py_body.append('residual[%i] = %s' % (ii, rhs)) 1699 c_rhs = ExprManip.make_c_compatible(rhs) 1700 c_body.append('residual[%i] = %s;' % (ii, c_rhs)) 1701 1702 py_body.append('') 1703 py_body.append('return residual') 1704 py_body = '\n '.join(py_body) 1705 1706 c_body.append('}') 1707 c_body = os.linesep.join(c_body) 1708 1709 self._dynamic_funcs_python.set('alg_res_func', py_body) 1710 self._dynamic_funcs_c.set('alg_res_func', c_body)

1711

1712 - def _make_dres_dc_function(self):

1713 py_body = [] 1714 py_body.append('def dres_dc_function(time, dynamicVars, yprime, ' 1715 'constants):') 1716 py_body.append('dynamicVars = scipy.asarray(dynamicVars)') 1717 py_body.append('yprime = scipy.asarray(yprime)') 1718 py_body.append('constants = scipy.asarray(constants)') 1719 py_body.append('') 1720 py_body.append('pd = scipy.zeros((%i, %i), scipy.float_)' 1721 %(len(self.dynamicVars), len(self.dynamicVars))) 1722 py_body.append('') 1723 self._add_assignments_to_function_body(py_body) 1724 py_body.append('') 1725 1726 c_body = [] 1727 c_args = 'double *time_ptr, double *dynamicVars, double *yprime, '\ 1728 'double *constants, double *pd' 1729 c_body.append('void dres_dc_function_(%s){' % c_args) 1730 self._prototypes_c.set('dres_dc_function', 1731 'void dres_dc_function_(%s);' % c_args) 1732 c_body.append('double time = *time_ptr;') 1733 c_body.append('') 1734 self._add_assignments_to_function_body(c_body, in_c=True) 1735 c_body.append('') 1736 1737 N_dyn = len(self.dynamicVars) 1738 1739 for res_ii, (rhs, dt) in enumerate(self._residual_terms): 1740 if rhs is None: 1741 rhs = '0' 1742 rhs_vars_used = ExprManip.extract_vars(rhs) 1743 for wrt_ii, wrt in enumerate(self.dynamicVars.keys()): 1744 pd_terms = [] 1745 deriv = self.takeDerivative(rhs, wrt, rhs_vars_used) 1746 if (deriv != '0' and deriv != '0.0'): 1747 py_body.append('# Residual %i wrt %s' % (res_ii, wrt)) 1748 py_body.append('pd[%i, %i] = %s' % (res_ii, wrt_ii, deriv)) 1749 deriv_c = ExprManip.make_c_compatible(deriv) 1750 c_body.append('pd[%i] = %s;' 1751 % (res_ii+wrt_ii*N_dyn, deriv_c)) 1752 1753 py_body.append('') 1754 py_body.append('return pd') 1755 py_body = '\n '.join(py_body) 1756 1757 c_body.append('}') 1758 c_body = os.linesep.join(c_body) 1759 1760 self._dynamic_funcs_python.set('dres_dc_function', py_body) 1761 self._dynamic_funcs_c.set('dres_dc_function', c_body)

1762

1763 - def _make_dres_dcdot_function(self):

1764 py_body = [] 1765 py_body.append('def dres_dcdot_function(time, dynamicVars, yprime, ' 1766 'constants):') 1767 py_body.append('dynamicVars = scipy.asarray(dynamicVars)') 1768 py_body.append('yprime = scipy.asarray(yprime)') 1769 py_body.append('constants = scipy.asarray(constants)') 1770 py_body.append('') 1771 py_body.append('pd = scipy.zeros((%i, %i), scipy.float_)' 1772 %(len(self.dynamicVars), len(self.dynamicVars))) 1773 py_body.append('') 1774 self._add_assignments_to_function_body(py_body) 1775 py_body.append('') 1776 1777 c_body = [] 1778 c_args = 'double *time_ptr, double *dynamicVars, double *yprime, '\ 1779 'double *constants, double *pd' 1780 c_body.append('void dres_dcdot_function_(%s){' % c_args) 1781 self._prototypes_c.set('dres_dcdot_function', 1782 'void dres_dcdot_function_(%s);' % c_args) 1783 c_body.append('double time = *time_ptr;') 1784 c_body.append('') 1785 self._add_assignments_to_function_body(c_body, in_c=True) 1786 c_body.append('') 1787 1788 N_dyn = len(self.dynamicVars) 1789 for res_ii, (rhs, dt) in enumerate(self._residual_terms): 1790 for wrt_ii, wrt in enumerate(self.dynamicVars.keys()): 1791 if dt == wrt: 1792 py_body.append('# Derivative of residual term %i wrt ' 1793 '%s_dot' % (res_ii, wrt)) 1794 py_body.append('pd[%i, %i] = -1' 1795 % (res_ii, wrt_ii)) 1796 c_body.append('pd[%i] = -1;' % (res_ii+wrt_ii*N_dyn)) 1797 1798 py_body.append('') 1799 py_body.append('return pd') 1800 py_body = '\n '.join(py_body) 1801 1802 c_body.append('}') 1803 c_body = os.linesep.join(c_body) 1804 1805 self._dynamic_funcs_python.set('dres_dcdot_function', py_body) 1806 self._dynamic_funcs_c.set('dres_dcdot_function', c_body)

1807

1808 - def _make_ddaskr_jac(self):

1809 py_body = [] 1810 py_body.append('def ddaskr_jac(time, dynamicVars, yprime, cj, ' 1811 'constants):') 1812 py_body.append('dynamicVars = scipy.asarray(dynamicVars)') 1813 py_body.append('yprime = scipy.asarray(yprime)') 1814 py_body.append('constants = scipy.asarray(constants)') 1815 py_body.append('') 1816 py_body.append('dres_dc = dres_dc_function(time, dynamicVars, yprime, ' 1817 'constants)') 1818 py_body.append('dres_dcdot = dres_dcdot_function(time, dynamicVars, ' 1819 'yprime, constants)') 1820 py_body.append('return dres_dc + cj * dres_dcdot') 1821 py_body = '\n '.join(py_body) 1822 1823 N_dyn = len(self.dynamicVars) 1824 c_body = [] 1825 c_args = 'double *time_ptr, double *dynamicVars, double *yprime, '\ 1826 'double *pd, double *cj_ptr, double *constants, int *intpar' 1827 c_body.append('void ddaskr_jac_(%s){' % c_args) 1828 self._prototypes_c.set('ddaskr_jac', 1829 'void ddaskr_jac_(%s);' % c_args) 1830 c_body.append('double cj = *cj_ptr;') 1831 c_body.append('') 1832 c_body.append('dres_dc_function_(time_ptr, dynamicVars, yprime, ' 1833 'constants, pd);') 1834 c_body.append('') 1835 # Like magic, this will initialize *all* elements to zero 1836 c_body.append('double dres_dcdot[%i*%i] = {0};' % (N_dyn, N_dyn)) 1837 c_body.append('dres_dcdot_function_(time_ptr, dynamicVars, yprime, ' 1838 'constants, dres_dcdot);') 1839 c_body.append('') 1840 c_body.append('int ii;') 1841 c_body.append('for(ii=0; ii < %i; ii++){' % N_dyn**2) 1842 c_body.append(' pd[ii] += cj*dres_dcdot[ii];}') 1843 c_body.append('}') 1844 c_body = os.linesep.join(c_body) 1845 1846 self._dynamic_funcs_python.set('ddaskr_jac', py_body) 1847 self._dynamic_funcs_c.set('ddaskr_jac', c_body)

1848

1849 - def _make_dres_dsinglep(self):

1850 # loop over the optimizable params, and for each one create 1851 # its own dres_dparam function 1852 for wrt_ii, wrt in enumerate(self.optimizableVars.keys()): 1853 1854 func_name = 'dres_d' + wrt 1855 1856 py_body = [] 1857 py_body.append('def ' + func_name + '(time, dynamicVars, ' 1858 'yprime, constants):') 1859 py_body.append('dynamicVars = scipy.asarray(dynamicVars)') 1860 py_body.append('constants = scipy.asarray(constants)') 1861 py_body.append('') 1862 py_body.append('pd = scipy.zeros(%i, scipy.float_)' 1863 % len(self.dynamicVars)) 1864 py_body.append('') 1865 1866 self._add_assignments_to_function_body(py_body) 1867 py_body.append('') 1868 1869 c_body = [] 1870 c_args = 'double *time_ptr, double *dynamicVars, double *yprime, '\ 1871 'double *constants, double *pd' 1872 c_body.append('void ' + func_name + '_(%s){' % c_args) 1873 self._prototypes_c.set(func_name, 1874 'void ' + func_name + '_(%s);' % c_args) 1875 1876 c_body.append('double time = *time_ptr;') 1877 c_body.append('') 1878 self._add_assignments_to_function_body(c_body, in_c=True) 1879 c_body.append('') 1880 1881 N_dyn = len(self.dynamicVars) 1882 N_opt = len(self.optimizableVars) 1883 1884 # Now the sensitivities. 1885 # We'll cache lists of the variables in each rhs, 1886 # since extract_vars is slow. 1887 rhs_vars = {} 1888 for (rhs, dt) in self._residual_terms: 1889 if rhs is None: 1890 rhs = '0' 1891 rhs_vars[rhs] = ExprManip.extract_vars(rhs) 1892 1893 for res_ii, (rhs, dt) in enumerate(self._residual_terms): 1894 if rhs is None: 1895 rhs = '0' 1896 rhs_vars_used = rhs_vars[rhs] 1897 deriv = self.takeDerivative(rhs, wrt, rhs_vars_used) 1898 if (deriv != '0' and deriv != '0.0'): 1899 py_body.append('pd[%i] = %s' % (res_ii, deriv)) 1900 c_deriv = ExprManip.make_c_compatible(deriv) 1901 c_body.append('pd[%i] = %s;' % (res_ii, c_deriv)) 1902 1903 py_body.append('') 1904 py_body.append('return pd') 1905 py_body = '\n '.join(py_body) 1906 1907 c_body.append('}') 1908 c_body = os.linesep.join(c_body) 1909 1910 self._dynamic_funcs_python.set(func_name, py_body) 1911 self._dynamic_funcs_c.set(func_name, c_body)

1912

1913 - def _make_sens_rhs(self):

1914 1915 py_body = [] 1916 py_body.append('def sens_rhs(time, sens_y, sens_yp, constants):') 1917 py_body.append('sens_y = scipy.asarray(sens_y)') 1918 py_body.append('sens_yp = scipy.asarray(sens_yp)') 1919 py_body.append('constants = scipy.asarray(constants)') 1920 py_body.append('') 1921 py_body.append('sens_res = scipy.zeros(%i, scipy.float_)' 1922 % (2*len(self.dynamicVars))) 1923 py_body.append('') 1924 1925 c_body = [] 1926 c_args = 'double *time_ptr, double *sens_y, double *sens_yp, '\ 1927 'double *cj_ptr, double *sens_res, int *ires_ptr, '\ 1928 'double *constants, int *ipar' 1929 c_body.append('void sens_rhs_(%s){' % c_args) 1930 self._prototypes_c.set('sens_rhs', 1931 'void sens_rhs_(%s);' % c_args) 1932 c_body.append('') 1933 1934 N_dyn = len(self.dynamicVars) 1935 N_const = len(self.constantVars) 1936 py_body.append('sens_res[:%(N_dyn)i] = res_function(time, ' 1937 'sens_y, sens_yp, constants)' 1938 % {'N_dyn': N_dyn}) 1939 py_body.append('') 1940 1941 # for the python sens_rhs function, we add a dictionary 1942 # of dres_dparam function names indexed by parameter index. 1943 py_body.append('func_dict = {}') 1944 for wrt_ii, wrt in enumerate(self.optimizableVars.keys()): 1945 py_body.append('func_dict[' + str(wrt_ii) + '] = dres_d' + wrt) 1946 1947 py_body.append('') 1948 1949 py_body.append('p_index = int(constants[%i])' % N_const) 1950 1951 # sens_rhs will access a particular dres_dparam function 1952 # based on the passed in parameter index 1953 py_body.append('dres_dp_func = func_dict[p_index]') 1954 1955 py_body.append('') 1956 1957 py_body.append('dc_dp = sens_y[%i:]' % N_dyn) 1958 py_body.append('dcdot_dp = sens_yp[%i:]' % N_dyn) 1959 py_body.append('dres_dp = dres_dp_func(time, sens_y, sens_yp, constants)' 1960 % {'N_dyn': N_dyn, 'N_const': N_const}) 1961 1962 py_body.append('dres_dc = dres_dc_function(time, ' 1963 'sens_y, sens_yp, constants)' 1964 % {'N_dyn': N_dyn, 'N_const': N_const}) 1965 py_body.append('dres_dcdot = dres_dcdot_function(time, ' 1966 'sens_y, sens_yp, constants)' 1967 % {'N_dyn': N_dyn, 'N_const': N_const}) 1968 py_body.append('sens_res[%i:] = dres_dp ' 1969 '+ scipy.dot(dres_dc, dc_dp) ' 1970 '+ scipy.dot(dres_dcdot, dcdot_dp)' % N_dyn) 1971 1972 # This fills in the first half of our sens_res 1973 c_body.append('res_function_(time_ptr, sens_y, sens_yp, cj_ptr, ' 1974 'sens_res, ires_ptr, constants, ipar);') 1975 c_body.append('') 1976 c_body.append('int p_index = (int)constants[%i];' % N_const) 1977 # we add an array that tracks the constants only, so that 1978 # the c version of the dres_dparam function gets an array of the 1979 # expected size 1980 c_body.append('double constants_only[%i];' % N_const) 1981 c_body.append('int jj;') 1982 c_body.append('for (jj = 0; jj < %s; jj++){' % N_const) 1983 c_body.append('constants_only[jj] = constants[jj];}') 1984 c_body.append('double *dc_dp = &sens_y[%i];' % N_dyn) 1985 c_body.append('double *dcdot_dp = &sens_yp[%i];' % N_dyn) 1986 # We'll directly fill dres_dp into the appropriate place in sens_res 1987 c_body.append('double *dres_dp = &sens_res[%i];' % N_dyn) 1988 c_body.append('int ii;') 1989 # sens_res isn't necessarily all zeros when passed in using the 1990 # cpointer. 1991 c_body.append('for(ii = 0; ii < %s; ii++){' % N_dyn) 1992 c_body.append('dres_dp[ii] = 0;}') 1993 1994 # we add a switch-case statement to choose the proper dres_dparam 1995 # function based on the passed in parameted index 1996 c_body.append('switch(p_index)') 1997 c_body.append('{') 1998 for wrt_ii, wrt in enumerate(self.optimizableVars.keys()): 1999 c_body.append('case ' + str(wrt_ii) + ' : dres_d' + wrt + '_(time_ptr, sens_y, sens_yp, constants_only, ' 2000 'dres_dp);') 2001 c_body.append('break;'), 2002 c_body.append('}') 2003 # Fill in dres_dc 2004 c_body.append('double dres_dc[%i] = {0};' % N_dyn**2) 2005 c_body.append('dres_dc_function_(time_ptr, sens_y, sens_yp, constants, ' 2006 'dres_dc);') 2007 c_body.append('int row, col;') 2008 c_body.append('for(row = 0; row < %i; row++){' % N_dyn) 2009 c_body.append('for(col = 0; col < %i; col++){' % N_dyn) 2010 c_body.append('sens_res[row+%i] += dres_dc[row + col*%i]*dc_dp[col];}}' 2011 % (N_dyn, N_dyn)) 2012 2013 c_body.append('double dres_dcdot[%i] = {0};' % N_dyn**2) 2014 c_body.append('dres_dcdot_function_(time_ptr, sens_y, sens_yp, ' 2015 'constants, dres_dcdot);') 2016 c_body.append('for(row = 0; row < %i; row++){' % N_dyn) 2017 c_body.append('for(col = 0; col < %i; col++){' % N_dyn) 2018 c_body.append('sens_res[row+%i] += dres_dcdot[row + col*%i]' 2019 '*dcdot_dp[col];}}' % (N_dyn, N_dyn)) 2020 2021 py_body.append('') 2022 py_body.append('return sens_res') 2023 py_body = '\n '.join(py_body) 2024 2025 c_body.append('}') 2026 c_body = os.linesep.join(c_body) 2027 2028 self._dynamic_funcs_python.set('sens_rhs', py_body) 2029 self._dynamic_funcs_c.set('sens_rhs', c_body)

2030

2031 - def _make_log_funcs(self):

2032 N_dyn = len(self.dynamicVars) 2033 2034 py_body = [] 2035 py_body.append('def res_function_logdv(time, log_dv, log_yp, ' 2036 'constants):') 2037 py_body.append('log_dv = scipy.asarray(log_dv)') 2038 py_body.append('log_yp = scipy.asarray(log_yp)') 2039 py_body.append('dynamicVars = scipy.exp(log_dv)') 2040 py_body.append('dynamicVars = scipy.maximum(dynamicVars, %g)' 2041 % _double_tiny_) 2042 py_body.append('yprime = log_yp * dynamicVars') 2043 py_body.append('return res_function(time, dynamicVars, yprime, ' 2044 'constants)') 2045 py_body = '\n '.join(py_body) 2046 self._dynamic_funcs_python.set('res_function_logdv', py_body) 2047 2048 c_body = [] 2049 c_args = 'double *time_ptr, double *log_dv, double *log_yp, '\ 2050 'double *cj_ptr, double *residual, int *ires_ptr, '\ 2051 'double *constants, int *ipar' 2052 c_body.append('void res_function_logdv_(%s){' % c_args) 2053 self._prototypes_c.set('res_function_logdv', 2054 'void res_function_logdv_(%s);' % c_args) 2055 c_body.append('double dynamicVars[%i];' % N_dyn) 2056 c_body.append('double yprime[%i];' % N_dyn) 2057 c_body.append('int ii;') 2058 c_body.append('for(ii = 0; ii < %i; ii++){' % N_dyn) 2059 c_body.append('dynamicVars[ii] = max(exp(log_dv[ii]), DBL_MIN);') 2060 c_body.append('yprime[ii] = log_yp[ii] * dynamicVars[ii];}') 2061 c_body.append('res_function_(time_ptr, dynamicVars, yprime, cj_ptr, ' 2062 'residual, ires_ptr, constants, ipar);') 2063 c_body.append('}') 2064 c_body = os.linesep.join(c_body) 2065 self._dynamic_funcs_c.set('res_function_logdv', c_body) 2066 2067 py_body = [] 2068 py_body.append('def root_func_logdv(time, log_dv, log_yp, constants):') 2069 py_body.append('log_dv = scipy.asarray(log_dv)') 2070 py_body.append('log_yp = scipy.asarray(log_yp)') 2071 py_body.append('dynamicVars = scipy.exp(log_dv)') 2072 py_body.append('dynamicVars = scipy.maximum(dynamicVars, %g)' 2073 % _double_tiny_) 2074 py_body.append('yprime = log_yp * dynamicVars') 2075 py_body.append('return root_func(time, dynamicVars, yprime, ' 2076 'constants)') 2077 py_body = '\n '.join(py_body) 2078 self._dynamic_funcs_python.set('root_func_logdv', py_body) 2079 2080 c_body = [] 2081 c_args = 'int *neq_ptr, double *time_ptr, double *log_dv, '\ 2082 'double *log_yp, int *nrt_ptr, double *root_devs, '\ 2083 'double *constants, int *ipar' 2084 c_body.append('void root_func_logdv_(%s){' % c_args) 2085 self._prototypes_c.set('root_func_logdv', 2086 'void root_func_logdv_(%s);' % c_args) 2087 c_body.append('double dynamicVars[%i];' % N_dyn) 2088 c_body.append('double yprime[%i];' % N_dyn) 2089 c_body.append('int ii;') 2090 c_body.append('for(ii = 0; ii < %i; ii++){' % N_dyn) 2091 c_body.append('dynamicVars[ii] = max(exp(log_dv[ii]), DBL_MIN);') 2092 c_body.append('yprime[ii] = log_yp[ii] * dynamicVars[ii];}') 2093 c_body.append('root_func_(neq_ptr, time_ptr, dynamicVars, yprime, ' 2094 'nrt_ptr, root_devs, constants, ipar);') 2095 c_body.append('}') 2096 c_body = os.linesep.join(c_body) 2097 self._dynamic_funcs_c.set('root_func_logdv', c_body) 2098 2099 2100 py_body = [] 2101 py_body.append('def sens_rhs_logdv(time, sens_y, sens_yp, constants):') 2102 # We need to copy sens_y and sens_yp here, since we aren't allowed to 2103 # change their values. 2104 py_body.append('sens_y = scipy.array(sens_y)') 2105 py_body.append('sens_yp = scipy.array(sens_yp)') 2106 py_body.append('sens_y[:%i] = scipy.exp(sens_y[:%i])' % (N_dyn, N_dyn)) 2107 py_body.append('sens_y[:%i] = scipy.maximum(sens_y[:%i], %g)' 2108 % (N_dyn, N_dyn, _double_tiny_)) 2109 py_body.append('sens_yp[:%i] = sens_yp[:%i] * sens_y[:%i]' 2110 % (N_dyn, N_dyn, N_dyn)) 2111 py_body.append('return sens_rhs(time, sens_y, sens_yp, constants)') 2112 py_body = '\n '.join(py_body) 2113 self._dynamic_funcs_python.set('sens_rhs_logdv', py_body) 2114 2115 c_body = [] 2116 c_args = 'double *time_ptr, double *sens_y_log, double *sens_yp_log, '\ 2117 'double *cj_ptr, double *sens_res, int *ires_ptr, '\ 2118 'double *constants, int *ipar' 2119 c_body.append('void sens_rhs_logdv_(%s){' % c_args) 2120 self._prototypes_c.set('sens_rhs_logdv_', 2121 'void sens_rhs_logdv_(%s);' % c_args) 2122 c_body.append('double sens_y[%i];' % (N_dyn*2)) 2123 c_body.append('double sens_yp[%i];' % (N_dyn*2)) 2124 c_body.append('int ii;') 2125 c_body.append('for(ii = 0; ii < %i; ii++){' % N_dyn) 2126 c_body.append('sens_y[ii] = max(exp(sens_y_log[ii]), DBL_MIN);') 2127 c_body.append('sens_yp[ii] = sens_yp_log[ii] * sens_y[ii];}') 2128 c_body.append('for(ii = %i; ii < %i; ii++){' % (N_dyn, 2*N_dyn)) 2129 c_body.append('sens_y[ii] = sens_y_log[ii];') 2130 c_body.append('sens_yp[ii] = sens_yp_log[ii];}') 2131 c_body.append('sens_rhs_(time_ptr, sens_y, sens_yp, cj_ptr, sens_res, ' 2132 'ires_ptr, constants, ipar);') 2133 c_body.append('}') 2134 c_body = os.linesep.join(c_body) 2135 self._dynamic_funcs_c.set('sens_rhs_logdv', c_body)

2136 2137

2138 - def _make_integrate_stochastic_tidbit(self):

2139 # Function definitions 2140 py_body = [] 2141 py_body.append('def integrate_stochastic_tidbit' 2142 '(seed, reseed, time, dv, cv, rmsd, stop_time):') 2143 2144 c_args = 'unsigned long* seed_ptr, int* reseed, double* time_ptr, '\ 2145 'int* dv, double* cv, '\ 2146 'double* rmsd_ptr, double* stop_time_ptr, double* trajectory' 2147 c_body = [] 2148 c_body.append('void integrate_stochastic_tidbit_(%s) {'%c_args) 2149 2150 # Test for conditions under which we cannot make these functions 2151 err_body = [] 2152 if len(self.reactions) == 0: 2153 err_body.append('# No reactions present.') 2154 if len(self.rateRules) > 0: 2155 err_body.append('# %i rate rules are present'%len(self.rateRules)) 2156 if len([_ for _ in self.events if not _.timeTriggered]) > 0: 2157 err_body.append('# %i non-time-triggered events are present'%\ 2158 len([_ for _ in self.events \ 2159 if not _.timeTriggered])) 2160 if len(self.algebraicRules) > 0: 2161 err_body.append('# %i algebraic rules are present'% 2162 len(self.algebraicRules)) 2163 if len(self.algebraicVars) > 0: 2164 err_body.append('# %i algebraic variables are present'% 2165 len(self.algebraicRules)) 2166 for rxn in self.reactions.values(): 2167 for val in rxn.stoichiometry.values(): 2168 try: intVal = int(val) 2169 except: intVal = None 2170 if intVal is None: 2171 err_body.append('# Noncastable stoichiometry for %s: %s'% 2172 (rxn.id, repr(val))) 2173 elif intVal != val: 2174 err_body.append('# Reaction stoichiometry for %s: %s!=%s'% 2175 (rxn.id, repr(val), repr(intVal))) 2176 2177 if err_body!=[]: 2178 py_body = py_body + err_body 2179 py_body.append('pass') 2180 py_body = '\n '.join(py_body) 2181 self._dynamic_funcs_python.set('integrate_stochastic_tidbit', 2182 py_body) 2183 2184 c_body.append('return;}') 2185 c_body = os.linesep.join(c_body) 2186 self._prototypes_c.set('integrate_stochastic_tidbit', 2187 'void integrate_stochatic_tidbit_(%s);' % 2188 c_args) 2189 self._dynamic_funcs_c.set('integrate_stochastic_tidbit', c_body) 2190 2191 return 2192 2193 N_CON = len(self.constantVars) 2194 N_DYN = len(self.dynamicVars) 2195 N_ASN = len(self.assignedVars) 2196 N_RXN = len(self.reactions) 2197 2198 # Use a custom AST walk to get 12x faster than ExprManip.sub_for_vars 2199 mapping = {} 2200 for index, name in enumerate(self.constantVars.keys()): 2201 mapping[name] = 'cv[%i]'%index 2202 for index, name in enumerate(self.dynamicVars.keys()): 2203 mapping[name] = 'dv[%i]'%index 2204 for index, name in enumerate(self.assignedVars.keys()): 2205 mapping[name] = 'av[%i]'%index 2206 2207 class Parse(ExprManip.AST.compiler.visitor.ASTVisitor): 2208 def __call__(slf,s,c=True): 2209 if c: s = ExprManip.make_c_compatible(s) 2210 ast = ExprManip.strip_parse(s) 2211 ExprManip.AST.compiler.walk(ast, slf) 2212 return ExprManip.ast2str(ast)

2213 2214 def visitName(slf,node,*args): 2215 if mapping.has_key(node.name): node.name=mapping[node.name] 2216 parse = Parse() 2217 2218 # Find the stoichiometry and reaction dependencies 2219 asnKeys = self.assignedVars.keys() 2220 evars = [] # The extracted dynamic variables each rxn depends upon 2221 for rxnInd, rxn in enumerate(self.reactions.values()): 2222 evar = ExprManip.Extraction.extract_vars(rxn.kineticLaw) 2223 evar = sets.Set(evar) 2224 while len(evar.intersection(asnKeys)): 2225 for asnName in evar.intersection(asnKeys): 2226 rule = self.assignmentRules.get(asnName) 2227 evar.remove(asnName) 2228 evar.union_update( \ 2229 ExprManip.Extraction.extract_vars(rule)) 2230 evars.append(evar) 2231 2232 stch_body = [] 2233 depd_body = [] 2234 for rxnInd, (rxn_id, rxn) in enumerate(self.reactions.items()): 2235 stch_body.append(repr([int(rxn.stoichiometry.get(_,0)) 2236 for _ in self.dynamicVars.keys()])[1:-1]) 2237 varNames = sets.Set([n 2238 for n,v in rxn.stoichiometry.items() if v!=0]) 2239 depd = [int(len(evar.intersection(varNames))>0) for evar in evars] 2240 depd_body.append(repr(depd)[1:-1]) 2241 if N_RXN>0: depd_body.append(repr([1]*N_RXN)[1:-1]) 2242 2243 # Python body assembly 2244 py_body.append('import numpy.random') 2245 py_body.append('if reseed: rs = numpy.random.seed(seed)') 2246 py_body.append('') 2247 if N_RXN > 0: 2248 py_body.append('stch = [[%s]]' % '],\n ['.join(stch_body)) 2249 py_body.append('depd = [[%s]]' % '],\n ['.join(depd_body)) 2250 else: 2251 py_body.append('stch = [[]]') 2252 py_body.append('depd = [[]]') 2253 py_body.append('') 2254 py_body.append('av = [0.0]*%i'%N_ASN) 2255 py_body.append('dv0 = [_ for _ in dv]') 2256 py_body.append('props = [0.0]*%i'%N_RXN) 2257 py_body.append('rxnInd = %i'%N_RXN) 2258 py_body.append('while time < stop_time:') 2259 for index, rule in enumerate(self.assignmentRules.values()): 2260 py_body.append(' av[%i]=%s;'%(index, parse(rule, False))) 2261 py_body.append('') 2262 for index, rxn in enumerate(self.reactions.values()): 2263 py_body.append(' if depd[rxnInd][%i]: props[%i]=%s'% 2264 (index, index, parse(rxn.kineticLaw, False))) 2265 py_body.append(' propensity=sum(props)') 2266 py_body.append(' if propensity==0.0:') 2267 py_body.append(' dt = stop_time - time') 2268 py_body.append(' time=stop_time') 2269 py_body.append(' break') 2270 py_body.append(' ') 2271 py_body.append(' dt = -log(1.0-numpy.random.rand())/propensity') 2272 py_body.append(' time += dt') 2273 py_body.append(' ') 2274 py_body.append(' selection=propensity*numpy.random.rand()') 2275 py_body.append(' for rxnInd in range(%i):'%N_RXN) 2276 py_body.append(' if selection<props[rxnInd]: break') 2277 py_body.append(' else: selection-=props[rxnInd]') 2278 py_body.append(' ') 2279 py_body.append(' for index, st in enumerate(stch[rxnInd]):') 2280 py_body.append(' dv[index] += st') 2281 py_body.append(' ') 2282 py_body.append(' _rmsd = sum([(_-__)**2 for _, __ in zip(dv,dv0)])') 2283 py_body.append(' _rmsd = (_rmsd/%i.)**0.5'%N_DYN) 2284 py_body.append(' if _rmsd>rmsd: break') 2285 py_body.append('') 2286 py_body.append('if time > stop_time: # Back interpolate') 2287 py_body.append(' trajectory = [' 2288 'float(_)+float(stch[rxnInd][ii])/dt*(stop_time-time) ' 2289 'for ii, _ in enumerate(dv)]') 2290 py_body.append('else:') 2291 py_body.append(' stop_time = time') 2292 py_body.append(' trajectory = [float(_) for _ in dv]') 2293 py_body.append('return time, dv, stop_time, trajectory') 2294 py_body = '\n '.join(py_body) 2295 2296 # C body assembly 2297 c_body.append(' int i; /* Temp variables */') 2298 c_body.append('') 2299 c_body.append(' unsigned long seed = *seed_ptr;') 2300 c_body.append(' if (*reseed) {init_genrand(seed);}') 2301 c_body.append('') 2302 if N_RXN > 0: 2303 c_body.append(' short stch[%i][%i] = {{%s}};' % 2304 (N_RXN, N_DYN, 2305 '},\n {'.join(stch_body))) 2306 c_body.append(' short depd[%i+1][%i] = {{%s}};' % 2307 (N_RXN, N_RXN, 2308 '},\n {'.join(depd_body))) 2309 else: 2310 c_body.append(' short stch[0][0];') 2311 c_body.append('') 2312 c_body.append(' double time = *time_ptr;') 2313 c_body.append(' double sd = (*rmsd_ptr)*(*rmsd_ptr)*%i.;'%N_DYN) 2314 c_body.append(' double stop_time = *stop_time_ptr;') 2315 c_body.append(' double dt=0.0;') 2316 c_body.append('') 2317 c_body.append(' double dv0[%i];'%N_DYN) 2318 c_body.append(' for (i=0;i<%i;i++) {dv0[i]=dv[i];}'%N_DYN) 2319 c_body.append('') 2320 c_body.append(' int rxnInd = %i;'%N_RXN) 2321 c_body.append(' double propensity, selection, props[%i], av[%i];'% 2322 (N_RXN, N_ASN)) 2323 c_body.append(' while (time < stop_time) {') 2324 for index, rule in enumerate(self.assignmentRules.values()): 2325 c_body.append(' av[%i]=%s;'%(index, parse(rule))) 2326 c_body.append('') 2327 for index, rxn in enumerate(self.reactions.values()): 2328 c_body.append(" if (depd[rxnInd][%i]) {props[%i]=%s;}"% \ 2329 (index, index, parse(rxn.kineticLaw))) 2330 c_body.append('') 2331 c_body.append(' propensity = 0.0;') 2332 c_body.append(' for (i=0;i<%i;i++) {'%N_RXN) 2333 c_body.append(' propensity += props[i];}') 2334 c_body.append(' if (propensity<=0.0) {') 2335 c_body.append(' dt = stop_time-time;') 2336 c_body.append(' time = stop_time;') 2337 c_body.append(' break;') 2338 c_body.append(' }') 2339 c_body.append('') 2340 c_body.append(' dt = -log(1.0-genrand_real32())/propensity;') 2341 c_body.append(' time += dt;') 2342 c_body.append('') 2343 c_body.append(' selection = propensity * genrand_real32();') 2344 c_body.append('') 2345 c_body.append(' for (rxnInd=0; rxnInd<%i; rxnInd++) {'%N_RXN) 2346 c_body.append(' if (selection < props[rxnInd]) {break;}') 2347 c_body.append(' else {selection -= props[rxnInd];}}') 2348 c_body.append('') 2349 c_body.append(' for (i=0;i<%i;i++) {dv[i]+=stch[rxnInd][i];}'%N_DYN) 2350 c_body.append('') 2351 c_body.append(' double _sd = 0.0;') 2352 c_body.append(' for (i=0;i<%i;i++) {'%N_DYN) 2353 c_body.append(' _sd += (dv0[i]-dv[i])*(dv0[i]-dv[i]);') 2354 c_body.append(' }') 2355 c_body.append(' if (_sd > sd) {break;}') 2356 c_body.append(' }') 2357 c_body.append('') 2358 c_body.append(' for (i=0;i<%i;i++) {'%N_DYN) 2359 c_body.append(' trajectory[i]=(double)dv[i];') 2360 c_body.append(' if (time > stop_time) {') 2361 c_body.append(' trajectory[i] += (double)stch[rxnInd][i]/dt*' 2362 '(stop_time - time);') 2363 c_body.append(' }') 2364 c_body.append(' }') 2365 c_body.append(' if (time>stop_time) {(*stop_time_ptr) = stop_time;}') 2366 c_body.append(' else {(*stop_time_ptr) = time;}') 2367 c_body.append('') 2368 c_body.append(' (*time_ptr) = time;') 2369 c_body.append('}') 2370 c_body.append('') 2371 c_body = os.linesep.join(c_body) 2372 2373 self._prototypes_c.set('integrate_stochastic_tidbit', 2374 'void integrate_stochatic_tidbit_(%s);' % c_args) 2375 2376 self._dynamic_funcs_python.set('integrate_stochastic_tidbit', py_body) 2377 self._dynamic_funcs_c.set('integrate_stochastic_tidbit', c_body) 2378 2379

2380 - def _add_assignments_to_function_body(self, body, in_c = False, 2381 include_dts=False):

2382 """ 2383 Adds the assignment rules for this Network to the list of lines that 2384 are in body. 2385 """ 2386 # We make sure these guys are arrays, so that we can be confident 2387 # in using .item(). 2388 # We loop to assign our constantVarValues and our dynamicVars 2389 # to local variables for speed (avoid repeated accesses) and 2390 # for readability 2391 for arg, var_names in zip(['constants', 'dynamicVars'], 2392 [self.constantVars.keys(), 2393 self.dynamicVars.keys()]): 2394 for ii, id in enumerate(var_names): 2395 if not in_c: 2396 body.append('%s = %s.item(%i)' % (id, arg, ii)) 2397 if include_dts and arg != 'constants': 2398 body.append('%s_deriv_wrt_time = yprime.item(%i)' 2399 % (id,ii)) 2400 else: 2401 body.append('double %s = %s[%i];' % (id, arg, ii)) 2402 if include_dts and arg != 'constants': 2403 body.append('double %s_deriv_wrt_time = yprime[%i];' 2404 % (id,ii)) 2405 body.append('') 2406 2407 for variable, math in self.assignmentRules.items(): 2408 if not in_c: 2409 body.append('%s = %s' % (variable, math)) 2410 else: 2411 c_math = ExprManip.make_c_compatible(math) 2412 body.append('double %s = %s;' % (variable, c_math)) 2413 if include_dts: 2414 # This zero is in case assigment doesn't depend on any 2415 # dynamic variables. In that case, simply_expr could fail 2416 # on an empty string. 2417 rhs_terms = ['0'] 2418 # Which variables does assignment rule depend on? 2419 dep_vars = ExprManip.extract_vars(math) 2420 2421 # Do the chain rule for dynamic variables. 2422 for dynvar in dep_vars.intersection(self.dynamicVars.keys()): 2423 dmath_dvar = ExprManip.diff_expr(math, dynvar) 2424 rhs_terms.append('%s * %s_deriv_wrt_time' 2425 % (dmath_dvar, dynvar)) 2426 # Add in explicit time dependence. 2427 if 'time' in dep_vars: 2428 dmath_dtime = ExprManip.diff_expr(math, 'time') 2429 rhs_terms.append(dmath_dtime) 2430 2431 rhs = ' + '.join(rhs_terms) 2432 rhs = ExprManip.simplify_expr(rhs) 2433 if not in_c: 2434 body.append('%s_deriv_wrt_time = %s' % (variable, rhs)) 2435 else: 2436 c_rhs = ExprManip.make_c_compatible(rhs) 2437 body.append('double %s_deriv_wrt_time = %s;' 2438 % (variable, c_rhs)) 2439 2440 return body

2441 2442 # 2443 # Methods involved in events 2444 #

2445 - def fireEvent(self, event, dynamicVarValues, time):

2446 self.updateVariablesFromDynamicVars(dynamicVarValues, time) 2447 delay = self.evaluate_expr(event.delay, time) 2448 2449 return delay

2450

2451 - def executeEvent(self, event, time_fired, y_fired, y_current, time_current):

2452 # The values assigned depend on the state of the network at the time 2453 # the event was >fired<. We go back to that time here... 2454 self.updateVariablesFromDynamicVars(y_fired, time_fired) 2455 new_values = {} 2456 var_vals = [(id, self.get_var_val(id)) for id in self.variables.keys()] 2457 var_vals = dict(var_vals) 2458 var_vals['time'] = time_fired 2459 # Calculate the values that will get assigned to the variables 2460 for id, rhs in event.event_assignments.items(): 2461 new_values[id] = self.evaluate_expr(rhs, time_fired, var_vals) 2462 2463 # Make all the relevant assignments 2464 for id, value in new_values.items(): 2465 y_current[self.dynamicVars.indexByKey(id)] = value 2466 # Update our network with the new values, just for consistency. 2467 self.updateVariablesFromDynamicVars(y_current, time_current) 2468 2469 return y_current

2470

2471 - def _smooth_trigger(self, trigger):

2472 ast = ExprManip.strip_parse(trigger) 2473 if len(ast.ops) != 1: 2474 raise ValueError('Comparison has more than one operation in ' 2475 'clause %s!' % trigger) 2476 lhs = ExprManip.ast2str(ast.expr) 2477 rhs = ExprManip.ast2str(ast.ops[0][1]) 2478 if ast.ops[0][0] in ['>', '>=']: 2479 return '%s - %s' % (lhs, rhs) 2480 elif ast.ops[0][0] in ['<', '<=']: 2481 return '-%s + %s' % (lhs, rhs) 2482 else: 2483 raise ValueError('Comparison %s in triggering clause is not ' 2484 '<, <=, >, or >=!')

2485 2486 # We cache derivatives wrt dynamicVariables, because they'll be redundantly 2487 # taken many, many times otherwise, and it can be a major hold-up. 2488 # Unforunately, the structures we want to index by aren't hashable, so this 2489 # can't be a dictionary. Luckily that shouldn't slow things down much, 2490 # because we shouldn't have to sort through many structures in our 2491 # applications. 2492 _sens_event_derivs_cache = []

2493 - def executeEventAndUpdateSens(self, holder, ysens_pre_exec, opt_var):

2494 curr_struct = self._get_structure() 2495 curr_events = copy.deepcopy(self.events) 2496 # Search through the _sens_event_derivs_cache for the derivs 2497 # corresponding to this structure. 2498 for struct, events, derivs_cache in self._sens_event_derivs_cache: 2499 if curr_struct == struct and curr_events == events: 2500 break 2501 else: 2502 derivs_cache = {} 2503 self._sens_event_derivs_cache.append((curr_struct, curr_events, 2504 derivs_cache)) 2505 2506 # We extract all our relevant quantities from the event_info holder 2507 event = holder.event 2508 time_fired = holder.time_fired 2509 ysens_fired = holder.ysens_fired 2510 yp_fired = holder.yp_fired 2511 clause_index = holder.clause_index 2512 time_exec = holder.time_exec 2513 yp_pre_exec = holder.yp_pre_exec 2514 y_post_exec = holder.y_post_exec 2515 yp_post_exec = holder.yp_post_exec 2516 2517 # This is subtle... If this event is the result of a chain, then the 2518 # firing time for the event is determined not by its trigger, but 2519 # by the execution time of the event that caused it to fire. 2520 # We walk up the list of chained events, so that finally we set the 2521 # d(time of firing) for this event to d(time of execution) for the 2522 # start of that chain of events. 2523 dtf_dp = None 2524 link = holder 2525 while link.chained_off_of is not None: 2526 dtf_dp = link.chained_off_of.dte_dp 2527 link = link.chained_off_of 2528 2529 ysens_post_exec = copy.copy(ysens_pre_exec) 2530 # Fill in the new values for the normal y variables 2531 N_dv = len(self.dynamicVars) 2532 ysens_post_exec[:N_dv] = y_post_exec 2533 2534 # We set our network to the state at which the event fired, so we 2535 # can use evaluate_expr for the calculations of d(firing time) 2536 self.updateVariablesFromDynamicVars(ysens_fired, time_fired) 2537 var_vals_fired = [(id, self.get_var_val(id)) for id 2538 in self.variables.keys()] 2539 var_vals_fired = dict(var_vals_fired) 2540 var_vals_fired['time'] = time 2541 # We need to do this if we don't have a chained event. 2542 if dtf_dp is None: 2543 # We convert our trigger to smooth function we can take the 2544 # derivative of so we can calculate sensitivity of firing time. 2545 # The trigger is passed to us because it might be a clause of a 2546 # more complicated trigger. 2547 trigger = self.event_clauses[clause_index] 2548 trigger = self._smooth_trigger(trigger) 2549 trigger_vars_used = ExprManip.extract_vars(trigger) 2550 2551 # Compute the sensitivity of our firing time 2552 # 2553 # numerator = sum_dv dtrigger_ddv ddv_dp + dtrigger_dp 2554 # denominator = sum_dv dtrigger_ddv ddv_dt + dtrigger_dt 2555 # dtf_dp = -(numerator)/(denominator) 2556 numerator = 0 2557 denominator = 0 2558 # First term in numerator: sum_dv dtrigger_ddv ddv_dp 2559 # First term in denominator: sum_dv dtrigger_ddv ddv_dt 2560 for ii, id in enumerate(self.dynamicVars.keys()): 2561 try: 2562 dtrigger_ddv = derivs_cache[(trigger, id)] 2563 except KeyError: 2564 dtrigger_ddv = self.takeDerivative(trigger, id, 2565 trigger_vars_used) 2566 derivs_cache[(trigger, id)] = dtrigger_ddv 2567 dtrigger_ddv = self.evaluate_expr(dtrigger_ddv, time_fired, 2568 var_vals_fired) 2569 # This is the index of ddv_dp in the y array 2570 index_in_y = ii + N_dv 2571 ddv_dp = ysens_fired[index_in_y] 2572 ddv_dt = yp_fired[ii] 2573 numerator += dtrigger_ddv * ddv_dp 2574 denominator += dtrigger_ddv * ddv_dt 2575 2576 # Second term in numerator: dtrigger_dp 2577 dtrigger_dp = self.takeDerivative(trigger, opt_var, 2578 trigger_vars_used) 2579 dtrigger_dp = self.evaluate_expr(dtrigger_dp, time_fired, 2580 var_vals_fired) 2581 numerator += dtrigger_dp 2582 2583 # Second term in denominator: dtrigger_dt 2584 dtrigger_dt = self.takeDerivative(trigger, 'time', 2585 trigger_vars_used) 2586 dtrigger_dt = self.evaluate_expr(dtrigger_dt, time_fired, 2587 var_vals_fired) 2588 denominator += dtrigger_dt 2589 2590 dtf_dp = -numerator/denominator 2591 2592 # Now compute the sensitivity of the delay 2593 # 2594 # ddelay/dp = ddelay/dp + ddelay_dy*dy/dp + ddelay_dy*dy/dt*dtf/dp 2595 # + ddelay/dt*dtf/dp 2596 delay = event.delay 2597 if not isinstance(delay, str): 2598 ddelay_dp = 0 2599 else: 2600 # First term 2601 delay_vars_used = ExprManip.extract_vars(delay) 2602 ddelay_dp = self.takeDerivative(delay, opt_var, delay_vars_used) 2603 ddelay_dp = self.evaluate_expr(ddelay_dp, time_fired, 2604 var_vals_fired) 2605 2606 # Second and third terms: ddelay_dy*dy/dp + ddelay_dy*dy/dt*dtf/dp 2607 for ii, id in enumerate(self.dynamicVars.keys()): 2608 try: 2609 ddelay_dy = derivs_cache[(delay, id)] 2610 except KeyError: 2611 ddelay_dy = self.takeDerivative(delay, id, delay_vars_used) 2612 derivs_cache[(delay, id)] = ddelay_dy 2613 ddelay_dy = self.evaluate_expr(ddelay_dy, time_fired, 2614 var_vals_fired) 2615 # This is the index of ddv_dp in the y array 2616 index_in_y = ii + N_dv 2617 dy_dp = ysens_fired[index_in_y] 2618 dy_dt = yp_fired[ii] 2619 ddelay_dp += ddelay_dy * dy_dp 2620 ddelay_dp += ddelay_dy * dy_dt * dtf_dp 2621 2622 # Fourth term: ddelay/dt*dtf/dp 2623 ddelay_dt = self.takeDerivative(delay, 'time', delay_vars_used) 2624 ddelay_dt = self.evaluate_expr(ddelay_dt, time_fired, 2625 var_vals_fired) 2626 ddelay_dp += ddelay_dt * dtf_dp 2627 2628 # This is the delay in the time of execution 2629 dte_dp = dtf_dp + ddelay_dp 2630 2631 # We store d(time of execution) for use by other events that may be 2632 # chained to this one. 2633 holder.dte_dp = dte_dp 2634 2635 # We update our Network so we can get all the variable values 2636 # that were just prior to the event executing. 2637 self.updateVariablesFromDynamicVars(ysens_pre_exec, time_exec) 2638 var_vals_pre_exec = [(id, self.get_var_val(id)) for id 2639 in self.variables.keys()] 2640 var_vals_pre_exec = dict(var_vals_pre_exec) 2641 var_vals_pre_exec['time'] = time 2642 2643 # Now compute the sensitivity of each of our new values 2644 for y_ii, y_id in enumerate(self.dynamicVars.keys()): 2645 if not event.event_assignments.has_key(y_id): 2646 # This is the index of y's sensitivity in the sensitivity array 2647 index_in_y = y_ii + N_dv 2648 # dy_dp after the event of course begins equal to the current 2649 # sensitivity 2650 dy_dp_post_exec = ysens_pre_exec[index_in_y] 2651 # These next two terms account for possible changes in the 2652 # derivative of y due to event execution. 2653 # dy/dt(pre_exec) * dte/dp term 2654 dy_dp_post_exec += yp_pre_exec[y_ii] * dte_dp 2655 # -dy/dt(post_exec) * dte/dp 2656 dy_dp_post_exec -= yp_post_exec[y_ii] * dte_dp 2657 ysens_post_exec[index_in_y] = dy_dp_post_exec 2658 else: 2659 # If y has been assigned the value a, then the sensitivity of 2660 # y after the event is da/dp - dy/dt(post_exec) * dte/dp 2661 a = event.event_assignments.get(y_id) 2662 if not isinstance(a, str): 2663 # If a isn't a string, it must be a constant, so da/dp = 0. 2664 dy_dp_post_exec = -yp_post_exec[y_ii] * dte_dp 2665 index_in_y = y_ii + N_dv 2666 ysens_post_exec[index_in_y] = dy_dp_post_exec 2667 else: 2668 # da/dp= da/dp + da/dy * dy/dp + da/dy * dy/dt * dtf/dp 2669 # + da/dt * dtf/dp 2670 # All these are calculated at the firing time 2671 2672 # We only need to deal with whatever clause in the 2673 # piecewise is applicable. 2674 a = self._sub_for_piecewise(a, time_fired) 2675 a_vars_used = ExprManip.extract_vars(a) 2676 2677 # First term: da/dp 2678 da_dp = self.takeDerivative(a, opt_var, a_vars_used) 2679 da_dp = self.evaluate_expr(da_dp, time_fired, 2680 var_vals_fired) 2681 dy_dp_post_exec = da_dp 2682 2683 # Second term: da/dy * dy/dp 2684 # Third term: da/dy * dy/dt * dtf/dp 2685 # Together = da/dy * (dy/dp + dy/dt * dtf/dp) 2686 for other_ii,other_id in enumerate(self.dynamicVars.keys()): 2687 try: 2688 da_dother = derivs_cache[(a, other_id)] 2689 except KeyError: 2690 da_dother = self.takeDerivative(a, other_id, 2691 a_vars_used) 2692 derivs_cache[(a, other_id)] = da_dother 2693 da_dother = self.evaluate_expr(da_dother, time_fired, 2694 var_vals_fired) 2695 2696 # This is the index of dother_dp in the y array 2697 index_in_y = other_ii + N_dv 2698 dother_dp = ysens_fired[index_in_y] 2699 2700 dother_dt = yp_fired[other_ii] 2701 2702 dy_dp_post_exec += da_dother * dother_dp 2703 dy_dp_post_exec += da_dother * dother_dt * dtf_dp 2704 2705 # Fourth term: da/dt * dtf/dp 2706 da_dt = self.takeDerivative(a, 'time', a_vars_used) 2707 da_dt = self.evaluate_expr(da_dt, time_fired, 2708 var_vals_fired) 2709 dy_dp_post_exec += da_dt * dtf_dp 2710 2711 # Finally, -dy/dt(post_exec) * dte/dp 2712 # We calculated dy/dt(new) up above in this function 2713 dy_dp_post_exec -= yp_post_exec[y_ii] * dte_dp 2714 2715 index_in_y = y_ii + N_dv 2716 ysens_post_exec[index_in_y] = dy_dp_post_exec 2717 2718 return ysens_post_exec

2719 2720 # 2721 # Internally useful things 2722 #

2723 - def _makeCrossReferences(self):

2724 """ 2725 Create the cross-reference lists for the Network. 2726 """ 2727 self.assignedVars = KeyedList() 2728 self.constantVars = KeyedList() 2729 self.optimizableVars = KeyedList() 2730 self.dynamicVars = KeyedList() 2731 self.algebraicVars = KeyedList() 2732 2733 self.compartments = KeyedList() 2734 self.parameters = KeyedList() 2735 self.species = KeyedList() 2736 mapping = {Compartment: self.compartments, 2737 Parameter: self.parameters, 2738 Species: self.species} 2739 2740 for id, var in self.variables.items(): 2741 mapping[var.__class__].set(id, var) 2742 if var.is_constant: 2743 self.constantVars.set(id, var) 2744 if var.is_optimizable: 2745 self.optimizableVars.set(id, var) 2746 elif id in self.assignmentRules.keys(): 2747 self.assignedVars.set(id, var) 2748 else: 2749 self.dynamicVars.set(id, var) 2750 2751 self.constantVarValues = [self.evaluate_expr(var.value) for var in 2752 self.constantVars.values()] 2753 self.constantVarValues = scipy.array(self.constantVarValues) 2754 2755 # Collect all variables that are explicitly in algebraic rules 2756 vars_in_alg_rules = sets.Set() 2757 for rule in self.algebraicRules: 2758 vars_in_alg_rules.union_update(ExprManip.extract_vars(rule)) 2759 2760 # Now replace all the assigned variables with the variables they 2761 # actually depend on. This takes a while loop because assigned 2762 # vars may be functions of other assigned vars. 2763 assigned_in_alg = vars_in_alg_rules.intersection(sets.Set( 2764 self.assignedVars.keys())) 2765 # while there are still assignment variables that we have not 2766 # expanded in terms of their definitions 2767 while assigned_in_alg: 2768 # Replace that assigned_var with the variables it depends on 2769 # in the running list of algebraic rules 2770 for assigned_var in assigned_in_alg: 2771 vars_in_alg_rules.union_update(ExprManip.extract_vars( 2772 self.assignmentRules.get(assigned_var))) 2773 vars_in_alg_rules.remove(assigned_var) 2774 # update the list of assignment variables that appear in the 2775 # algebraic rule 2776 assigned_in_alg = vars_in_alg_rules.intersection( 2777 sets.Set(self.assignedVars.keys())) 2778 2779 # At this point, vars_in_alg_rules should contain all the variables the 2780 # algebraic rules depend on, including implicit dependencies through 2781 # assignment rules. Now we filter out all the things we know aren't 2782 # algebraic vars. First we filter out everything that we already 2783 # know isn't a dynamic variable. 2784 vars_in_alg_rules.intersection_update(sets.Set(self.dynamicVars.keys())) 2785 2786 # remove the reaction variables 2787 for rxn in self.reactions: 2788 for chem, value in rxn.stoichiometry.items(): 2789 if value != 0 and (chem in vars_in_alg_rules): 2790 vars_in_alg_rules.remove(chem) 2791 2792 # remove the rate variables 2793 for var in self.rateRules.keys(): 2794 if (var in vars_in_alg_rules): 2795 vars_in_alg_rules.remove(var) 2796 2797 # remove the event variables 2798 for e in self.events: 2799 for var in e.event_assignments.keys(): 2800 if (var in vars_in_alg_rules): 2801 vars_in_alg_rules.remove(var) 2802 2803 # Set the algebraicVars list to the list we just compiled, after sorting 2804 # based on the order in self.dynamicVars 2805 sorted_alg_vars = [(id, self.variables.get(id)) 2806 for id in self.dynamicVars.keys() 2807 if id in vars_in_alg_rules] 2808 self.algebraicVars = KeyedList(sorted_alg_vars)

2809 2810

2811 - def disable_deriv_funcs(self):

2812 """ 2813 Disables the creation of derivative functions, which can speed up 2814 integration. 2815 """ 2816 2817 if self.deriv_funcs_enabled == True: 2818 self.ddaskr_jac = None 2819 self._dynamic_structure_methods = ['_make_res_function', 2820 '_make_alg_deriv_func', 2821 '_make_alg_res_func', 2822 '_make_integrate_stochastic_tidbit' 2823 ] 2824 self.deriv_funcs_enabled = False

2825

2826 - def enable_deriv_funcs(self):

2827 """ 2828 Enables the creation of derivative functions, which allows 2829 sensitivity integration and may speed up normal integration 2830 as well. 2831 """ 2832 2833 if self.deriv_funcs_enabled == False: 2834 self._dynamic_structure_methods = ['_make_res_function', 2835 '_make_alg_deriv_func', 2836 '_make_alg_res_func', 2837 '_make_dres_dc_function', 2838 '_make_dres_dcdot_function', 2839 '_make_ddaskr_jac', 2840 '_make_dres_dsinglep', 2841 '_make_sens_rhs', 2842 '_make_log_funcs', 2843 '_make_integrate_stochastic_tidbit' 2844 ] 2845 self.deriv_funcs_enabled = True

2846 2847 2848 _last_structure = None 2849 # This is an option to disable compilation of C modules. 2850 disable_c = SloppyCell.disable_c 2851 _last_disabled_c = disable_c

2852 - def compile(self, disable_c=None):

2853 """ 2854 Create the dynamically-generated functions for this Network. 2855 2856 Note that if a model involves many copies of the same network, with 2857 differences only in events or initial conditions, it will save time 2858 to compile the base Network before copying it. 2859 """ 2860 # If the structure of our network has changed, remake all the dynamic 2861 # functions 2862 # Note that this runs at least once, since _last_structure doesn't 2863 # exist beforehand. Also note that __setstate__ will exec all our 2864 # dynamic functions upon copying. 2865 if disable_c is None: 2866 disable_c = self.disable_c 2867 2868 curr_structure = self._get_structure() 2869 structure_changed = (curr_structure != self._last_structure) 2870 if structure_changed: 2871 self._last_structure = copy.deepcopy(curr_structure) 2872 # The dynamic function lists need to be reset because the number 2873 # and names of the dres_dparam functions can differ between 2874 # networks. 2875 # We do not include the list of function definitions because 2876 # they are handled differently than other dynamic functions. 2877 # Users should make sure that if they remove a function definition 2878 # from a compiled network that they also remove references to 2879 # that function from the rest of the network 2880 2881 # We clear out all the dynamic functions that have been defined. 2882 all_dynamic_keys = sets.Set(self._dynamic_funcs_python.keys()) 2883 all_dynamic_keys.union_update(self._dynamic_funcs_c.keys()) 2884 for dynamic_func in all_dynamic_keys: 2885 try: 2886 delattr(self, dynamic_func) 2887 except AttributeError: 2888 pass 2889 2890 self._dynamic_funcs_python = KeyedList() 2891 self._prototypes_c = self._common_prototypes_c.copy() 2892 self._dynamic_funcs_c = KeyedList() 2893 2894 reexec = False 2895 if structure_changed: 2896 logger.debug('Network %s: compiling structure.' % self.id) 2897 self._makeCrossReferences() 2898 # Check whether system appears well-determined. 2899 if len(self.algebraicVars) > len(self.algebraicRules): 2900 raise ValueError('System appears under-determined. ' 2901 'Not enough algebraic rules. ' 2902 '\n algebraicVars (%i) = \n%s ' 2903 '\n algebraicRules (%i) = \n%s '%(len(self.algebraicVars), self.algebraicVars, 2904 len(self.algebraicRules), self.algebraicRules)) 2905 self.make_func_defs() 2906 self._makeDiffEqRHS() 2907 for method in self._dynamic_structure_methods: 2908 getattr(self, method)() 2909 reexec = True 2910 2911 if self.events != getattr(self, '_last_events', None)\ 2912 or structure_changed: 2913 logger.debug('Network %s: compiling events.' % self.id) 2914 self._makeCrossReferences() 2915 for method in self._dynamic_event_methods: 2916 getattr(self, method)() 2917 self._last_events = copy.deepcopy(self.events) 2918 reexec = True 2919 2920 # after compile, check that there are not more algebraic variables than 2921 # algebraic rules. 2922 if len(self.algebraicVars) > len(self.algebraicRules): 2923 raise ValueError('System appears under-determined. ' 2924 'Not enough algebraic rules. ' 2925 '\n algebraicVars (%i) = \n%s ' 2926 '\n algebraicRules (%i) = \n%s '%(len(self.algebraicVars), self.algebraicVars, 2927 len(self.algebraicRules), self.algebraicRules)) 2928 2929 if self._last_disabled_c != disable_c: 2930 self._last_disabled_c = disable_c 2931 reexec = True 2932 2933 self.compiled=True 2934 if reexec: 2935 self.exec_dynamic_functions(disable_c=disable_c)

2936

2937 - def _get_structure(self):

2938 """ 2939 Return a tuple representing the structure of the Network. 2940 2941 The tuple contains the functionDefinitions, reactions, and rules for 2942 the Network. It also contains information on the constancy and 2943 optimizability of the variables. 2944 """ 2945 var_struct = {} 2946 structure = (self.functionDefinitions, self.reactions, 2947 self.assignmentRules, self.rateRules, var_struct, 2948 self.algebraicRules) 2949 for id, var in self.variables.items(): 2950 # If a constant variable is set equal to a function of other 2951 # variables, we should include that function, otherwise 2952 # our sensitivities will be wrong. 2953 if isinstance(self.get_var_ic(id), str): 2954 var_struct[id] = (var.is_constant, var.is_optimizable, 2955 self.get_var_ic(id)) 2956 else: 2957 var_struct[id] = (var.is_constant, var.is_optimizable) 2958 2959 return structure

2960 2961 # We cache our exec'd python functions and compiled C modules to minimize 2962 # redundant compiles. 2963 _py_func_dict_cache = {} 2964 _c_module_cache = {} 2965 # This is an option to disable compilation of C modules.

2966 - def exec_dynamic_functions(self, disable_c=False, del_c_files=True, 2967 curr_c_code=None):

2968 # only get the bodies that were created. 2969 curr_py_bodies = '\n'.join([body for body in self._dynamic_funcs_python.values()\ 2970 if body != None]) 2971 2972 key = (curr_py_bodies, tuple(self._func_strs.items())) 2973 # Search our cache of python functions. 2974 try: 2975 py_func_dict = self._py_func_dict_cache[key] 2976 except KeyError: 2977 # We don't have a cached version, so we need to generate it. 2978 py_func_dict = {} 2979 self._py_func_dict_cache[key] = py_func_dict 2980 for func_name, body in self._dynamic_funcs_python.items(): 2981 if body != None: 2982 exec body in self.namespace, locals() 2983 py_func_dict[func_name] = locals()[func_name] 2984 2985 # Add all the functions to our Network. 2986 for func_name, func in py_func_dict.items(): 2987 setattr(self, func_name, func) 2988 self.namespace[func_name] = func 2989 2990 if disable_c: 2991 return 2992 2993 if curr_c_code is None: 2994 curr_c_code = self.get_c_code() 2995 # Search the cache. 2996 try: 2997 c_module = self._c_module_cache[curr_c_code] 2998 except KeyError: 2999 # Regenerate if needed. 3000 # Write C to file. 3001 module_name = self.output_c(curr_c_code) 3002 try: 3003 # Run f2py on the C. This may raise an exception if the command 3004 # fails. 3005 self.run_f2py(module_name, hide_f2py_output=True) 3006 c_module = __import__(module_name) 3007 if del_c_files: 3008 os.unlink('%s.pyf' % module_name) 3009 os.unlink('%s.c' % module_name) 3010 try: 3011 os.unlink('%s.so' % module_name) 3012 except OSError: 3013 pass 3014 try: 3015 os.unlink('%s.pyd' % module_name) 3016 except OSError: 3017 pass 3018 except ImportError, X: 3019 # Compiling C failed. 3020 logger.warn('Failed to import dynamically compiled C module %s ' 3021 'for network %s.' % (module_name, self.get_id())) 3022 logger.warn(X) 3023 # We stored None for the c_module, so we don't repeatedly 3024 # try compiling the same bad C code. 3025 c_module = None 3026 self._c_module_cache[curr_c_code] = c_module 3027 3028 # Now we add all the appropriate functions to our Network. 3029 if c_module is not None: 3030 self.import_c_funcs_from_module(c_module)

3031

3032 - def get_c_code(self):

3033 # Combine all our C functions into one block of code 3034 c_code = [] 3035 c_code.append('#include <math.h>') 3036 c_code.append('#include <stdio.h>') 3037 c_code.append('#include <float.h>') 3038 c_code.append('#include "mtrand.h"') 3039 c_code.append('#define exponentiale M_E') 3040 c_code.append('#define pi M_PI') 3041 c_code.append('double max(double a, double b){') 3042 c_code.append('return a > b ? a : b;}') 3043 3044 # Function prototypes 3045 for func_name, proto in self._prototypes_c.items(): 3046 c_code.append(proto) 3047 c_code.append('') 3048 # Functions necessary for SBML math support 3049 for function, body in self._common_func_strs_c.items(): 3050 c_code.append(body) 3051 c_code.append('') 3052 # Function definitions 3053 for func_name, body in self._func_defs_c.items(): 3054 c_code.append(body) 3055 c_code.append('') 3056 # The dynamic functions for this network 3057 for func_name, body in self._dynamic_funcs_c.items(): 3058 c_code.append(body) 3059 c_code.append('') 3060 3061 c_code = '\n'.join(code for code in c_code if code != None) 3062 return c_code

3063

3064 - def output_c(self, c_code, mod_name=None):

3065 logger.debug('Outputting C for network %s.' % self.get_id()) 3066 if mod_name is None: 3067 semi_unique = str(time.time()) 3068 mod_name = '%s_%i_%s' % (self.get_id(), SloppyCell.my_rank, 3069 semi_unique[::-1]) 3070 mod_name = mod_name.replace('-', '_') 3071 # f2py fails if base filename has a '.' in it. 3072 mod_name = mod_name.replace('.', '_') 3073 3074 # Write the C code to a file. 3075 c_fd = open('%s.c' % mod_name, 'w') 3076 c_fd.write(c_code) 3077 c_fd.close() 3078 3079 # Read in the signatures that we'll fill in 3080 # use different sig file if derivs are disabled 3081 if self.deriv_funcs_enabled != True: 3082 pyf_base_filename = os.path.join(SloppyCell.__path__[0], 3083 'ReactionNetworks', 3084 'f2py_signatures_no_derivs.pyf') 3085 else: 3086 pyf_base_filename = os.path.join(SloppyCell.__path__[0], 3087 'ReactionNetworks', 3088 'f2py_signatures.pyf') 3089 3090 pyf_base_fd = file(pyf_base_filename, 'r') 3091 pyf_base = pyf_base_fd.read() 3092 3093 pyf_base_fd.close() 3094 3095 N_dyn = len(self.dynamicVars) 3096 N_const = len(self.constantVars) 3097 N_alg = len(self.algebraicVars) 3098 3099 # Fill in the signatures. 3100 # We begin by generating the signatures for the dres_dparam code. 3101 # Since the number of optimizable parameters varies in each network, 3102 # we cannot include these functions in the f2p_signatres template. 3103 dres_dparams_code = '' 3104 if self.deriv_funcs_enabled == True: 3105 for wrt_ii, wrt in enumerate(self.optimizableVars.keys()): 3106 func_name = 'dres_d' + wrt 3107 dres_dparams_code += ' subroutine ' + func_name + '(time, dynamicVars, yprime, constants, pd)\n' 3108 dres_dparams_code += ' double precision intent(in) :: time\n' 3109 dres_dparams_code += ' double precision intent(in), dimension(' + str(N_dyn) + ') :: dynamicVars\n' 3110 dres_dparams_code += ' double precision intent(in), dimension(' + str(N_dyn) + ') :: yprime\n' 3111 dres_dparams_code += ' double precision intent(in), dimension(' + str(N_const) + ') :: constants\n' 3112 dres_dparams_code += ' double precision intent(out), dimension(' + str(N_dyn) + ') :: pd\n' 3113 dres_dparams_code += ' end subroutine ' + func_name + '\n' 3114 3115 # Now update the template with the code we just generated and the other 3116 # names and dimensions 3117 pyf_code = pyf_base % {'dres_dparams': dres_dparams_code, 3118 'mod_name': mod_name, 3119 'N_dyn': N_dyn, 3120 'N_const': N_const, 3121 'N_alg': N_alg, 3122 'N_rt': self.len_root_func} 3123 3124 # Write out the signature file. 3125 pyf_fd = open('%s.pyf' % mod_name, 'w') 3126 pyf_fd.write(pyf_code) 3127 pyf_fd.close() 3128 3129 return mod_name

3130

3131 - def run_f2py(self, mod_name, hide_f2py_output=True):

3132 logger.debug('Running f2py for network %s.' % self.get_id()) 3133 from numpy.distutils.exec_command import exec_command 3134 # Run f2py. The 'try... finally...' structure ensures that we stop 3135 # redirecting output if there's an exection in running f2py 3136 # These options assume we're working with mingw. 3137 win_options = '' 3138 if sys.platform == 'win32': 3139 win_options = '--compiler=mingw32 --fcompiler=gnu' 3140 try: 3141 if hide_f2py_output: 3142 redir = Utility.Redirector() 3143 redir.start() 3144 sc_path = os.path.join(SloppyCell.__path__[0], 'ReactionNetworks') 3145 command = '-c %(win_options)s %(mod_name)s.pyf '\ 3146 '%(mod_name)s.c %(sc_path)s/mtrand.c -I%(sc_path)s'\ 3147 % {'win_options': win_options, 'mod_name': mod_name, 3148 'sc_path': sc_path} 3149 3150 prev_cflags = os.getenv('CFLAGS', '') 3151 os.putenv('CFLAGS', prev_cflags + ' -Wno-unused-variable') 3152 # f2py wants an extra argument at the front here. It's not actually 3153 # used though... 3154 oldargv = sys.argv 3155 sys.argv = ['f2py'] + command.split() 3156 import numpy.f2py.f2py2e 3157 output = numpy.f2py.f2py2e.run_compile() 3158 sys.argv = sys.argv 3159 os.putenv('CFLAGS', prev_cflags) 3160 except SystemExit, X: 3161 logger.warn('Call to f2py failed for network %s.' % self.get_id()) 3162 logger.warn(X) 3163 finally: 3164 if hide_f2py_output: 3165 redir.stop()

3166

3167 - def import_c_funcs_from_module(self, module):

3168 for function in self._dynamic_funcs_c.keys(): 3169 setattr(self, function, getattr(module, function))

3170

3171 - def takeDerivative(self, input, wrt, vars_used=None, simplify=True):

3172 """ 3173 Take the derivative of a math expression wrt a given variable id. 3174 3175 Does the chain rule through assigned variables. 3176 """ 3177 output = ExprManip.diff_expr(input, wrt) 3178 3179 if vars_used is None: 3180 vars_used = ExprManip.extract_vars(input) 3181 3182 # What other assigned variables does input depend on? 3183 assigned_used = vars_used.difference(sets.Set([wrt])) 3184 assigned_used.intersection_update(sets.Set(self.assignedVars.keys())) 3185 # Do the chain rule for those variables 3186 for id in assigned_used: 3187 rule = self.assignmentRules.getByKey(id) 3188 d2 = self.takeDerivative(rule, wrt, simplify=False) 3189 if d2 != '0': 3190 d = ExprManip.diff_expr(input, id) 3191 output += ' + (%s) *(%s)' % (d, d2) 3192 3193 # What other constant variables does input depend on? 3194 constant_used = vars_used.difference(sets.Set([wrt])) 3195 constant_used.intersection_update(sets.Set(self.constantVars.keys())) 3196 # Do the chain rule for those variables 3197 for id in constant_used: 3198 ic = self.get_var_ic(id) 3199 if isinstance(ic, str): 3200 d2 = self.takeDerivative(ic, wrt, simplify=False) 3201 if d2 != '0': 3202 d = ExprManip.diff_expr(input, id) 3203 output += ' + (%s) *(%s)' % (d, d2) 3204 3205 if simplify and output != '0': 3206 return ExprManip.simplify_expr(output) 3207 else: 3208 return output

3209

3210 - def copy(self, new_id=None, new_name=None):

3211 """ 3212 Return a copy of the given network, with an optional new id. 3213 """ 3214 new_net = copy.deepcopy(self) 3215 if new_id is not None: 3216 new_net.set_id(new_id) 3217 if new_name is not None: 3218 new_net.set_name(new_name) 3219 3220 return new_net

3221

3222 - def __getstate__(self):

3223 # deepcopy automatically does a deepcopy of whatever we return 3224 # here, so we only need to do a shallow copy and remove functions 3225 odict = copy.copy(self.__dict__) 3226 3227 # We can't copy the functions themselves. So we strip them out. 3228 for func in self._dynamic_funcs_python.keys(): 3229 odict[func] = None 3230 for func in self._dynamic_funcs_c.keys(): 3231 odict[func] = None 3232 odict['namespace'] = None 3233 # Let's not pickle these since they can be large and it would slow 3234 # down parallel execution. 3235 odict['trajectory'] = None 3236 odict['ddv_dpTrajectory'] = None 3237 3238 return odict

3239

3240 - def __setstate__(self, newdict):

3241 self.__dict__.update(newdict) 3242 self.namespace = copy.copy(self._common_namespace) 3243 3244 # Recreate our namespace 3245 for func_id, func_str in self._func_strs.items(): 3246 self.namespace[func_id] = eval(func_str, self.namespace, {}) 3247 3248 self._makeCrossReferences() 3249 if self.compiled: 3250 self.exec_dynamic_functions(self._last_disabled_c)

3251

3252 - def get_component_name(self, id, TeX_form=False):

3253 """ 3254 Return a components's name if it exists, else just return its id. 3255 """ 3256 # These are all the things that have names (except the network itself) 3257 complists = [self.variables, self.reactions, self.functionDefinitions, 3258 self.events, self.constraints] 3259 # If we don't find a name, we'll just use the id 3260 name = id 3261 for complist in complists: 3262 # If the id is in a list and has a non-empty name 3263 if complist.has_key(id) and complist.get(id).name: 3264 name = complist.get(id).name 3265 break 3266 3267 # We can also check the network's name 3268 if id == self.id and self.name: 3269 name = self.name 3270 3271 if TeX_form: 3272 # If we've got one underscore in the name, use that to indicate a 3273 # subscript 3274 if name.count('_') == 1: 3275 sp = name.split('_') 3276 name = '%s_{%s}' % (sp[0], sp[1]) 3277 else: 3278 # TeX can't handle more than one _ in a name, so we substitute 3279 name = name.replace('_', r'\_') 3280 3281 return name

3282

3283 - def get_eqn_structure(self):

3284 # This was used to interface with PyDSTool. 3285 out = {} 3286 out['odes'] = dict(self.diff_eq_rhs.items()) 3287 out['functions'] = {} 3288 for func_id, func_def in self.functionDefinitions.items(): 3289 vars = ', '.join(func_def.variables) 3290 out['functions']['%s(%s)' % (func_id, vars)] = func_def.math 3291 out['parameters'] = dict([(id, var.value) for (id, var) 3292 in self.constantVars.items()]) 3293 out['assignments'] = dict(self.assignmentRules.items()) 3294 out['events'] = dict([(event.trigger, event.event_assignments) 3295 for event in self.events]) 3296 out['constraints'] = dict([constraint.trigger for constraint in self.constraints]) 3297 3298 return out

3299 3300 3301 # Deprecated functions below.

3302 - def addCompartment(self, id, size=1.0, name='', 3303 typicalValue=False,isConstant=True, isOptimizable=False):

3304 self.add_compartment(id = id, initial_size = size, name = name, 3305 typical_value = typicalValue, 3306 is_constant = isConstant, 3307 is_optimizable = isOptimizable)

3308 3309 addVariable = _add_variable 3310

3311 - def addSpecies(self, id, compartment, initialConcentration=None, 3312 name='', typicalValue=None, is_boundary_condition=False, 3313 isConstant=False, isOptimizable=False, uniprot_ids=set()):

3314 self.add_species(id = id, compartment = compartment, 3315 initial_conc = initialConcentration, 3316 name = name, 3317 typical_value = typicalValue, 3318 is_boundary_condition = is_boundary_condition, 3319 is_constant = isConstant, 3320 is_optimizable = isOptimizable, 3321 uniprot_ids=uniprot_ids)

3322

3323 - def addParameter(self, id, value = 0.0, 3324 typicalValue = None, name = '', 3325 isConstant = True, isOptimizable = True):

3326 self.add_parameter(id = id, initial_value = value, name = name, 3327 typical_value = typicalValue, 3328 is_constant = isConstant, 3329 is_optimizable = isOptimizable)

3330 3331 addRateRule = add_rate_rule 3332

3333 - def dyn_var_fixed_point(self, dv0 = None):

3334 return Dynamics.dyn_var_fixed_point(self, dv0)

3335 FindFixedPoint = dyn_var_fixed_point 3336 3337 set_initial_var_value = set_var_ic 3338 setInitialVariableValue = set_var_ic 3339 3340 setOptimizables = update_optimizable_vars 3341

3342 - def addEvent(self, id, trigger, eventAssignments, delay=0, name=''):

3343 self.add_event(id = id, trigger = trigger, 3344 event_assignments = eventAssignments, name = name, 3345 delay = delay)

3346

3347 - def addConstraint(self, id, trigger, message=None, name=''):

3348 self.add_constraint(id = id, trigger = trigger, message = message, 3349 name = name)

3350 3351 3352 addFunctionDefinition = add_func_def 3353 addAssignmentRule = add_assignment_rule 3354 3355 3356

3357 -def _exec_dynamic_func(obj, func, in_namespace={}, bind=True):

3358 """ 3359 Create the executable function corresponding to func's functionBody. 3360 3361 This only exists now for Trajectory_mod. It's not used by the Network 3362 object. 3363 """ 3364 try: 3365 function_body = obj._dynamic_funcs_python.get(func) 3366 except (KeyError, AttributeError): 3367 function_body = getattr(obj, '%s_functionBody' % func) 3368 # This exec gives the function access to everything defined in in_namespace 3369 # and inserts the result into the locals namespace 3370 exec function_body in in_namespace, locals() 3371 # The call to types.MethodType ensures that we can call the function 3372 # as obj.f(...) and get the implicit 'self' argument. 3373 # locals()[func] just gets the actual function object the exec created. 3374 # Note that this this does depend on the _functionBody using a def 3375 # with the proper name. 3376 setattr(obj, func, 3377 types.MethodType(locals()[func], obj, obj.__class__))

3378

Source Code for Module SloppyCell.ReactionNetworks.Network_mod