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