Source Code for Module libTI2.solvers.csolver

  1  #!/usr/bin/env python 
  2   
  3  import copy,popen2,os,time,datetime,subprocess 
  4   
  5  scipy2c = {"power":"pow"} 
  6   
  7 -class csolver(object): 
  8      "docstring for csolver" 
  9 -    def __init__(self, species, species_values, differential_equations, signalling, precision=1e-10,stepsize=0.00001): 
 10          """If the signalling flag is set, simulation will be performed twice to get the AUC and other signal characteristica.""" 
 11          self.signalling = signalling 
 12          self.species = species 
 13          self.species_values = species_values 
 14          self.differential_equations = differential_equations 
 15          self.precision=precision 
 16          self.stepsize = stepsize 
 17          self.equationstring = "" 
 18          self.initialvaluestring = "" 
 19          infile = file("csolver.c","r") 
 20          self.content = "".join( infile.readlines() ) 
 21          infile.close() 
 22          # prepare differential equations 
 23          for species_index in range(len(self.species)): 
 24              species = self.species[species_index] 
 25              for var in self.differential_equations: 
 26                  self.differential_equations[var] = self.differential_equations[var].replace(" "+species+" "," y["+str(species_index)+"] ") 
 27          for species_index in range(len(self.species)): 
 28              species = self.species[species_index] 
 29              self.equationstring += "  out["+str(species_index)+"] = "+self.differential_equations[species]+";\n" 
 30              self.initialvaluestring += "  y["+str(species_index)+"] = "+str(self.species_values[species])+";\n" 
 31          self.cleanequationstring() 
 32          self.replacies = {} 
 33          self.replacies["signalling_flag"] = "" if self.signalling else "//" 
 34          self.replacies["putyourequationshere"] = self.equationstring 
 35          self.replacies["putyourinitialdefinitionshere"] = self.initialvaluestring 
 36          self.replacies["problemsize"] = str(len(self.species)) 
 37          self.replacies["stepsize"] = str(stepsize) 
 38          self.replacies["precision"] = str(self.precision) 
 39 -    def cleanequationstring(self): 
 40          for s,c in scipy2c.iteritems(): 
 41              self.equationstring = self.equationstring.replace("scipy."+s,c) 
 42 -    def filetostring(self): 
 43          for i,o in self.replacies.iteritems(): 
 44              self.content = self.content.replace("/*"+i+"*/",o) 
 45          return self.content 
 46 -    def tofile(self,filename): 
 47          f = file(filename,"w") 
 48          f.write(self.filetostring()) 
 49          f.close() 
 50 -    def solve(self,filename=None): 
 51   
 52          self.steady_state = {} 
 53          self.tau = {} 
 54          self.sigma = {} 
 55          self.ampl = {} 
 56          self.auc = {} 
 57           
 58   
 59          if not filename: 
 60              filename = self.random_filename() 
 61          self.tofile(filename+".c") 
 62          command1 = "gcc "+filename+".c -o "+filename+" -lm" 
 63          command2 = "./"+filename 
 64   
 65          try: 
 66              (o,i)=popen2.popen2(command1) 
 67              for line in o: print "Error during compilation:",o 
 68              # found this code snippet on code.pui.ch 
 69              # execute code with a timeout set 
 70              timeout_in_secs = 200 
 71              beginning = datetime.datetime.now() 
 72              p = subprocess.Popen(command2.split(" "), stdout=subprocess.PIPE, stderr=subprocess.PIPE ) 
 73              while p.poll() == None: 
 74                  time.sleep(0.1) 
 75                  current = datetime.datetime.now() 
 76                  if (current-beginning).seconds > timeout_in_secs: 
 77                      # timeout of 200 seconds 
 78                      raise RuntimeError,"Timeout for the c solver" 
 79              o = p.stdout.readlines() 
 80   
 81              # parsing the output of the c program 
 82              for line in o: 
 83                  #print line.strip("\n") 
 84                  if not line.startswith("result"): continue 
 85                  entries = line.strip("\n").split(",")[1:] 
 86                  species = self.species[ int(entries[0]) ] 
 87                  self.steady_state[species] = float(entries[1]) 
 88                  self.tau[species] = float(entries[2]) 
 89                  self.sigma[species] = float(entries[3]) 
 90                  self.ampl[species] = float(entries[4]) 
 91                  self.auc[species] = float(entries[5]) 
 92          except RuntimeError: 
 93              import scipy,signal 
 94              os.kill( p.pid, signal.SIGKILL ) 
 95              os.waitpid( -1, os.WNOHANG ) 
 96              for s in self.species: 
 97                  for l in [self.steady_state,self.tau,self.sigma,self.ampl,self.auc]: 
 98                      l[s] = scipy.nan 
 99              print "Timeout on reaction system" 
100          finally: 
101              os.remove(filename+".c") 
102              try: 
103                  os.remove(filename) 
104              except IOError: 
105                  print "Could not remove file",filename 
106 -    def random_filename(self): 
107          import random,string 
108          return "TIde"+"".join( random.Random().sample(string.letters+string.digits,12) ) 
109   
110           
111   
112  if __name__ == '__main__': 
113      import datetime 
114      s = ["x","y","z"] 
115      ev = {"x":" 0.1 - 0.01 * x ","y":" x - 0.1 * y * y ","z":" x * x + y * y * x - 11 * z * z "} 
116      iv = {"x":1000,"y":100,"z":0} 
117      #s = ["x"] 
118      #ev = {"x":" - x "} 
119      #iv = {"x":1} 
120   
121      start = datetime.datetime.now() 
122      cs = csolver(s,iv,ev) 
123      cs.solve() 
124      print datetime.datetime.now() - start 
125      print "Steady state" 
126      print cs.steady_state 
127      print "Tau" 
128      print cs.tau 
129      print "Sigma" 
130      print cs.sigma 
131      print "Amplitude" 
132      print cs.ampl 
133      print "AUC" 
134      print cs.auc 
135