1 from __future__ import division
2
3 import logging
4 logger = logging.getLogger('RxnNets.Trajectory_mod')
5
6 import os
7 import copy
8 import types
9
10 import scipy
11 import scipy.io
12 import scipy.interpolate
13
14 import SloppyCell.KeyedList_mod
15 KeyedList = SloppyCell.KeyedList_mod.KeyedList
16 import Network_mod
17 import SloppyCell.ExprManip as ExprManip
18
20 _known_structures = []
21 _known_function_bodies = []
22 _dynamic_funcs = ['_assignment', '_sens_assignment']
23
24 _common_namespace = {'log': scipy.log,
25 'log10': scipy.log10,
26 'exp': scipy.exp,
27 'cos': scipy.cos,
28 'sin': scipy.sin,
29 'tan': scipy.tan,
30 'acos': scipy.arccos,
31 'asin': scipy.arcsin,
32 'atan': scipy.arctan,
33 'cosh': scipy.cosh,
34 'sinh': scipy.sinh,
35 'tanh': scipy.tanh,
36 'arccosh': scipy.arccosh,
37 'arcsinh': scipy.arcsinh,
38 'arctanh': scipy.arctanh,
39 'pow': scipy.power,
40 'sqrt': scipy.sqrt,
41 'exponentiale': scipy.e,
42 'pi': scipy.pi,
43 }
44
45 - def __init__(self, net, key_column=None, is_sens=False, holds_dt=False,
46 empty=False, const_vals=None):
47 if empty:
48 return
49
50 if key_column is not None:
51 self.key_column = key_column
52 else:
53 keys = net.dynamicVars.keys() + net.assignedVars.keys()
54 if is_sens:
55 keys.extend([(cname, pname) for cname in keys
56 for pname in net.optimizableVars.keys()])
57 if holds_dt:
58 for keyname in copy.copy(keys):
59 if isinstance(keyname,str):
60 keys.append((keyname,'time'))
61 else:
62 keys.append(keyname + ('time',))
63
64 self.key_column = KeyedList(zip(keys, range(len(keys))))
65
66
67 self.timepoints = scipy.zeros(0, scipy.float_)
68 self.values = scipy.zeros((0, len(self.key_column)), scipy.float_)
69
70 self.var_keys = net.variables.keys()
71 self.dynamicVarKeys = net.dynamicVars.keys()
72 self.assignedVarKeys = net.assignedVars.keys()
73 self.optimizableVarKeys = net.optimizableVars.keys()
74
75
76
77 if const_vals is None:
78 self.const_var_values = KeyedList([(id, net.evaluate_expr(id)) for
79 id in net.constantVars.keys()])
80 else:
81 self.const_var_values = KeyedList(zip(net.constantVars.keys(),
82 const_vals))
83
84 self.typical_var_values = KeyedList([(id, var.typicalValue)
85 for (id, var)
86 in net.variables.items()])
87 self.event_info = None
88 self.tcks = {}
89 self.dytcks = {}
90
91
92 self._func_strs = copy.copy(net._func_strs)
93 self.namespace = copy.copy(self._common_namespace)
94 for func_id, func_str in self._func_strs.items():
95 self.namespace[func_id] = eval(func_str, self.namespace, {})
96
97
98
99 curr_structure = (net._last_structure, self.key_column, is_sens)
100 for ii, struct in enumerate(self._known_structures):
101 if curr_structure == struct:
102 (self._assignment_functionBody,
103 self._sens_assignment_functionBody) = \
104 self._known_function_bodies[ii]
105 break
106 else:
107
108
109
110 self._assignment_functionBody = self._make__assignment(net)
111 if is_sens:
112 self._sens_assignment_functionBody =\
113 self._make__sens_assignment(net)
114
115 self._known_structures.append(copy.deepcopy(curr_structure))
116 bodies = (self._assignment_functionBody,
117 getattr(self, '_sens_assignment_functionBody', None))
118 self._known_function_bodies.append(bodies)
119
121 return len(self.timepoints)
122
124
125
126 new_traj = copy.deepcopy(self)
127 if not isinstance(this_slice, slice):
128
129 this_slice = slice(this_slice, this_slice+1, None)
130 new_traj.timepoints = self.timepoints[this_slice]
131 new_traj.values = self.values[this_slice]
132
133
134
135 if self.tcks or new_traj.dytcks:
136 logger.warn('Interpolating functions must be recreated after '
137 'slicing a trajectory. Could be fixed.')
138 new_traj.tcks = {}
139 new_traj.dytcks = {}
140
141 return new_traj
142
144 """
145 Return a new trajectory containing only times from start to stop.
146 """
147 start_index = self._get_time_index(start, eps)
148 stop_index = self._get_time_index(stop, eps)
149 return self[start_index:stop_index+1]
150
153
155 return self.timepoints
156
157 - def add_event_info(self, net, eventinfo, time_traj_ended, include_extra_event_info = False):
158 """
159 Add information about the network state at each event execution
160 """
161
162
163
164
165
166 if include_extra_event_info == True:
167 (te, ye_pre, ye_post, ie) = eventinfo
168 assigned_states = []
169 prev_vals = net.get_var_vals()
170 for t, y, ii in zip(te, ye_pre, ie):
171 net.updateVariablesFromDynamicVars(y, t)
172 a_ids = [id for id in net.assignedVars.keys()]
173 a_vals = [net.get_var_val(id) for id in net.assignedVars.keys()]
174
175 a_state = dict(zip(a_ids,a_vals))
176 assigned_states.append(a_state)
177 net.set_var_vals(prev_vals, time_traj_ended)
178 eventinfo = (te,ye_pre,ye_post,ie,assigned_states)
179
180
181
182 self.event_info = eventinfo
183
184
185
187 if self.key_column != other.key_column:
188 raise ValueError, 'Trajectories in append have different column keys!'
189 if self.const_var_values != other.const_var_values:
190 logger.warn('Constant variable values differ between appended trajectories!')
191
192 if self.timepoints[-1] > other.timepoints[0]:
193 logger.warn('Appending trajectory with earlier timepoints!')
194
195 self.timepoints = scipy.concatenate((self.timepoints, other.timepoints))
196 self.values = scipy.concatenate((self.values, other.values))
197
199 return self.typical_var_values.get(id)
200
202 if self.key_column.has_key(id):
203 return self.values[:, self.key_column.get(id)]
204 elif self.const_var_values.has_key(id):
205 return scipy.ones(len(self.timepoints), scipy.float_) *\
206 self.const_var_values.get(id)
207 elif id == 'time':
208 return self.get_times()
209 else:
210 raise ValueError, 'Variable %s not found in trajectory.' % str(id)
211
213 """
214 Return the index of the stored value closest to the requested time.
215
216 Prints a warning if the difference between the requested time and the
217 stored time is greater than a fraction eps of the trajectory length.
218 """
219 index = scipy.argmin(abs(self.timepoints - time))
220 time_range = self.timepoints[-1] - self.timepoints[0]
221 if abs(self.timepoints[index] - time)/time_range > eps:
222 logger.warn('Time %f requested, closest time stored in trajectory '
223 'is %f.' % (time, self.timepoints[index]))
224 return index
225
227 """
228 Return a KeyedList of the values of the trajectory's variables at the
229 given time.
230
231 Prints a warning if the difference between the requested time and the
232 stored time is greater than a fraction eps of the trajectory length.
233 """
234 index = self._get_time_index(time, eps)
235 return self.get_var_vals_index(index)
236
238 """
239 Return the value of the given variable at the given time.
240
241 Prints a warning if the difference between the requested time and the
242 stored time is greater than a fraction eps of the trajectory length.
243 """
244 index = self._get_time_index(time, eps)
245 return self.get_var_val_index(var_id, index)
246
248 """
249 Return a KeyedList of the values of the trajectory's variables at the
250 given index.
251 """
252 out = KeyedList([(key, self.get_var_val_index(key, index)) for
253 key in ['time'] + self.keys()])
254 return out
255
257 """
258 Return the value of the given variable at the given index.
259 """
260 if self.key_column.has_key(var_id):
261 col = self.key_column.get(var_id)
262 return self.values[index, col]
263 elif self.const_var_values.has_key(var_id):
264 return self.const_var_values.get(var_id)
265 elif var_id == 'time':
266 return self.timepoints[index]
267
269 """
270 Return a KeyedList of the values of the trajectory's dynamic variables
271 at the given index.
272 """
273 out = KeyedList([(key, self.get_var_val_index(key, index)) for
274 key in self.dynamicVarKeys])
275 return out
276
278 """
279 Return a KeyedList of the values of the trajectory's dynamic variables
280 at the given time.
281
282 Prints a warning if the difference between the requested time and the
283 stored time is greater than a fraction eps of the trajectory length.
284 """
285 index = self._get_time_index(time, eps)
286 return self.get_dynvar_vals_index(index)
287
289 functionBody = ['def _assignment(self, values, times, start, end):']
290
291 for id in self.const_var_values.keys():
292 functionBody.append("%s = self.const_var_values.get('%s')" %
293 (id, id))
294
295 if len(net.assignmentRules) > 0:
296 for id, rule in net.assignmentRules.items():
297 lhs = self._sub_var_names(id)
298
299 rhs = self._sub_var_names(rule)
300 functionBody.append('# Assignment rule %s = %s' % (id, rule))
301 functionBody.append('%s = %s' % (lhs, rhs))
302 else:
303 functionBody.append('pass')
304
305 return '\n\t'.join(functionBody) + '\n'
306
308 functionBody = ['def _sens_assignment(self, values, times, start, end):'
309 ]
310
311 for id in self.const_var_values.keys():
312 functionBody.append("%s = self.const_var_values.get('%s')" %
313 (id, id))
314
315 if len(net.assignmentRules) > 0:
316 for id, rule in net.assignmentRules.items():
317
318 derivWRTdv = {}
319 for wrtId in net.dynamicVars.keys():
320 deriv = net.takeDerivative(rule, wrtId)
321 if deriv != '0':
322 derivWRTdv[wrtId] = deriv
323
324 for optId in net.optimizableVars.keys():
325 lhs = self._sub_var_names('%s__derivWRT__%s' % (id,optId))
326 rhs = []
327
328
329 for wrtId, deriv in derivWRTdv.items():
330 rhs.append('(%s) * %s__derivWRT__%s' %
331 (deriv, wrtId, optId))
332
333
334 derivWRTp = net.takeDerivative(rule, optId)
335 if derivWRTp != '0':
336 rhs.append(derivWRTp)
337
338 if rhs:
339 rhs = ' + '.join(rhs)
340 rhs = self._sub_var_names(rhs)
341 functionBody.append('%s = %s' % (lhs, rhs))
342 else:
343 functionBody.append('pass')
344
345 return '\n\t'.join(functionBody) + '\n'
346
348 if getattr(self, '_assignment', None) is None:
349 Network_mod._exec_dynamic_func(self, '_assignment',
350 self.namespace, bind=False)
351
352 numAdded = odeint_array.shape[0]
353 addedValues = scipy.zeros((numAdded, len(self.key_column)),
354 scipy.float_)
355
356 self.values = scipy.concatenate((self.values, addedValues))
357 self.timepoints = scipy.concatenate((self.timepoints, timepoints))
358
359 for ii, id in enumerate(self.dynamicVarKeys):
360 self.values[-numAdded:, self.key_column.get(id)] =\
361 odeint_array[:, ii]
362
363 self._assignment(self.values, self.timepoints, -numAdded, None)
364 if holds_dt :
365 for ii, id in enumerate(self.dynamicVarKeys) :
366 self.values[-numAdded:, self.key_column.get((id,'time'))] = \
367 odeint_array[:,ii+len(self.dynamicVarKeys)]
368
370 if getattr(self, '_assignment', None) is None:
371 Network_mod._exec_dynamic_func(self, '_assignment',
372 self.namespace, bind=False)
373
374 if getattr(self, '_sens_assignment', None) is None:
375 Network_mod._exec_dynamic_func(self, '_sens_assignment',
376 self.namespace, bind=False)
377
378 numAdded = odeint_array.shape[0]
379 addedValues = scipy.zeros((numAdded, len(self.key_column)),
380 scipy.float_)
381
382 self.values = scipy.concatenate((self.values, addedValues))
383 self.timepoints = scipy.concatenate((self.timepoints, timepoints))
384
385 nDv = len(self.dynamicVarKeys)
386 nOv = len(self.optimizableVarKeys)
387
388
389 for ii, dvId in enumerate(self.dynamicVarKeys):
390 self.values[-numAdded:, self.key_column.get(dvId)] = \
391 odeint_array[:, ii]
392
393 for ii, dvId in enumerate(self.dynamicVarKeys):
394 for jj, ovId in enumerate(self.optimizableVarKeys):
395 self.values[-numAdded:,
396 self.key_column.get((dvId, ovId))]\
397 = odeint_array[:, ii + (jj+1)*nDv]
398
399 self._assignment(self.values, self.timepoints, -numAdded, None)
400 self._sens_assignment(self.values, self.timepoints, -numAdded, None)
401
402 if holds_dt :
403
404 for ii, dvId in enumerate(self.dynamicVarKeys):
405 self.values[-numAdded:, self.key_column.get((dvId,'time'))] = \
406 odeint_array[:, ii + nDv*(nOv+1)]
407
408 for ii, dvId in enumerate(self.dynamicVarKeys):
409 for jj, ovId in enumerate(self.optimizableVarKeys):
410 self.values[-numAdded:,
411 self.key_column.get((dvId, ovId,'time'))]\
412 = odeint_array[:, ii + (jj+1)*nDv + nDv*(nOv+1)]
413
415 """
416 Return a copy of this trajectory containing only the variables specified
417 in keys.
418
419 If keys is None, all variables are included.
420 """
421 if keys == None:
422 keys = self.key_column.keys()
423
424 state = self.__getstate__()
425
426
427
428 keys = [key for key in keys if self.key_column.has_key(key)]
429 new_key_column = KeyedList(zip(keys, range(len(keys))))
430 state['key_column'] = new_key_column
431
432 new_values = scipy.zeros((len(self.values), len(new_key_column)),
433 scipy.float_)
434 for key, new_col in new_key_column.items():
435 old_col = self.values[:, self.key_column.get(key)]
436 new_values[:, new_col] = old_col.copy()
437 state['values'] = new_values
438
439 new_traj = Trajectory(None, empty=True)
440 new_traj.__setstate__(state)
441 return new_traj
442
444 odict = copy.copy(self.__dict__)
445 odict['_assignment'] = None
446 odict['_sens_assignment'] = None
447 odict['namespace'] = None
448 return odict
449
451 self.__dict__.update(newdict)
452
453 self.namespace = copy.copy(self._common_namespace)
454 for func_id, func_str in self._func_strs.items():
455 self.namespace[func_id] = eval(func_str, self.namespace, {})
456
458 mapping_dict = {}
459 for id in ExprManip.extract_vars(input):
460
461
462
463
464 splitId = id.split('__derivWRT__')
465 if len(splitId) == 1:
466 idname = splitId[0]
467 elif len(splitId) == 2:
468 idname = tuple(splitId)
469 else:
470 raise 'Problem with id %s in Trajectory._sub_var_names' % id
471
472 if idname in self.key_column.keys():
473 mapping = 'values[start:end, %i]' % self.key_column.get(idname)
474 elif idname in self.const_var_values.keys():
475
476
477 continue
478 elif idname == 'time':
479 mapping = 'times[start:end]'
480 else:
481 raise 'Problem with idname %s in Trajectory._sub_var_names' % id
482 mapping_dict[id] = mapping
483
484 input = ExprManip.sub_for_vars(input, mapping_dict)
485
486 return input
487
489 """ Given that a trajectory exists, build_interpolated_traj will create
490 the coefficients for the spline interpolatation.
491 The spline can then be evaluated using
492 Trajectory.evaluate_interpolated_traj or
493 Trajectory.evaluate_interpolated_trajs """
494 te,ye,ie = self.event_info[:3]
495 teIndices = []
496
497 if len(te) == 0 :
498 intervals = [(0,len(self.timepoints))]
499 else :
500
501
502 last_t = -1
503 for tevent in te :
504 if tevent != last_t:
505 teIndices.append(scipy.nonzero(self.timepoints==tevent)[0][1])
506 last_t = tevent
507
508
509
510 teIndicesWith0 = list(teIndices[0:])
511 teIndicesWith0.insert(0,0)
512
513
514 teIndices.extend([len(self.timepoints)])
515 intervals = zip(teIndicesWith0,teIndices)
516
517 self.tcks = {}
518
519 for (start_ind,end_ind) in intervals :
520 start_time, end_time = self.timepoints[start_ind], \
521 self.timepoints[end_ind-1]
522 curTimes = self.timepoints[start_ind:end_ind]
523 k = min(5,end_ind-start_ind-1)
524 ys = [self.get_var_traj(dv_id)[start_ind:end_ind]
525 for dv_id in self.key_column.keys()]
526
527 self.tcks[(start_time,end_time)] = [scipy.interpolate.splrep(curTimes,scipy.asarray(y),k=k,s=0) for y in ys]
528
529
530
532 """
533 This is a version of evaluate_interpolated_traj that returns all the
534 values of the dynamic variables and requires you to pass in the
535 appropriate subinterval between events (that can be found in
536 Trajectory.tcks.keys() )
537 Faster than calling evaluate_interpolated_traj repeatedly
538 """
539 local_tcks = self.tcks
540
541 nDVs = len(self.dynamicVarKeys)
542 dv_y = [scipy.interpolate.splev(time, local_tcks[subinterval][dv_ind],
543 der=der) for dv_ind in range(0,nDVs)]
544
545 return dv_y
546
548 """ Needs Trajectory.build_interpolated_traj() to be called first
549
550 Arguments:
551 dvid the name of the component of the trajectory you wish to
552 evaluate
553 time a vector of times or a scalar
554 subinterval an optional argument specifying the time interval
555 between events that the time argument lies (but given a
556 valid time, it will be found automatically)
557 der the derivative of the spline function you want, the order
558 of the derivative will be constrained by the order of the
559 interpolated spline
560 Outputs:
561 A single scalar value (if time input is a scalar)
562 or
563 (returned_times, interpolated_trajectory at those times) if times is a
564 vector
565
566 Note: It is necessary to have a returned_times argument too, in case
567 the times passed in happens to have a timepoint that corresponds
568 to an event time, which often has two trajectory values associated
569 with it.
570 """
571 if scipy.isscalar(time) :
572 time = scipy.asarray([time])
573 else :
574 time = scipy.asarray(time)
575 local_tcks = self.tcks
576 sorted_intervals = scipy.sort(local_tcks.keys(),axis=0)
577
578 if subinterval is not None :
579 if subinterval not in local_tcks.keys() :
580 raise "Not a valid subinterval (not in Trajectory.tcks.keys())"
581 else :
582 sorted_intervals = [[subinterval[0],subinterval[1]]]
583 interval_start_ind = 0
584 interval_end_ind = 0
585 else :
586
587 for interval_ind, interval in enumerate(sorted_intervals) :
588 start_time, end_time = interval[0],interval[1]
589 if (time[0] >= start_time) :
590 interval_start_ind = interval_ind
591 if (time[-1] <= end_time) :
592 interval_end_ind = interval_ind
593 break
594
595 dv_y = []
596 returned_times = []
597 dv_ind = self.key_column.keyToIndex[dv_id]
598 for interval in sorted_intervals[interval_start_ind:(interval_end_ind+1)] :
599 currTimes = scipy.compress( scipy.logical_and((time>=interval[0]),(time<=interval[1])) , time )
600 startslice, endslice = 0, None
601 if len(currTimes) > 1 :
602 if (currTimes[0]==currTimes[1]) :
603
604 startslice = 1
605 if (currTimes[-1]==currTimes[-2]) :
606
607 endslice = -1
608 dv_y.extend( scipy.interpolate.splev(currTimes[startslice:endslice],
609 local_tcks[(interval[0],interval[1])][dv_ind],der=der) )
610 returned_times.extend(currTimes[startslice:endslice])
611 elif len(currTimes) == 1:
612 dv_y.extend( [ scipy.interpolate.splev(currTimes, local_tcks[(interval[0],interval[1])][dv_ind],der=der) ])
613 returned_times.extend(currTimes[startslice:endslice])
614
615 if len(returned_times) == 1 :
616 return dv_y[0]
617 else :
618 return returned_times,dv_y
619
620 - def to_file(self, file_name, out_vars=None, separator=', '):
621 """
622 Output the given variables to a file.
623
624 file_name Name of the file to use for output (will be overwritten!)
625 out_vars List of variable ids ot output. If None, default is
626 'time', dynamic variables
627 separator The separator to use between columns.
628 """
629 if out_vars is None:
630 out_vars = ['time'] + self.dynamicVarKeys
631
632 first_line = separator.join(out_vars) + os.linesep
633 f = file(file_name, 'w')
634 f.write(first_line)
635
636 out_array = []
637 for var in out_vars:
638 out_array.append(self.get_var_traj(var))
639
640 out_array = scipy.transpose(out_array)
641 f = scipy.io.write_array(f, out_array, separator=separator,
642 keep_open=True)
643 f.close()
644
645
646
648 """
649 Return a list of the dynamic variable values at the last timepoint in
650 the trajectory.
651 """
652 return [self.values[-1, self.key_column.get(dv_id)] for dv_id in
653 self.dynamicVarKeys]
654
655 getVariableTrajectory = get_var_traj
656
657 try:
658 import psyco
659 psyco.bind(scipy.interpolate.splrep)
660 psyco.bind(scipy.interpolate.splev)
661 psyco.bind(evaluate_interpolated_trajs)
662 psyco.bind(evaluate_interpolated_traj)
663 except ImportError:
664 pass
665