Source Code for Module libTI2.optimize.nnca

  1  #!/usr/bin/env python 
  2  import scipy,random,scipy.linalg,scipy.optimize,datetime,copy,pickle,math,pickle,numpy 
  3   
  4 -def fmin_brute(f,ranges,logarithmic=True): 
  5      if type(ranges) == numpy.ndarray: 
  6          newranges = [] 
  7          numpoints = 2 
  8          for i in range(len(ranges)): 
  9              newranges.append( [ranges[i]/2.0,ranges[i]*2.0,numpoints] ) 
 10          ranges = newranges 
 11      y_to_x = {} 
 12      def range_to_list(tuple): 
 13          if not logarithmic: 
 14              return [((tuple[1]-tuple[0])/tuple[2]*y+tuple[0]) for y in range(tuple[2]+1)] 
 15          else: 
 16              return [math.exp((math.log(tuple[1])-math.log(tuple[0]))/tuple[2]*y+math.log(tuple[0])) for y in range(tuple[2]+1)] 
 17      def fix_part(f,ranges,fixed_x): 
 18          if len(ranges) == 0: 
 19              y = f(scipy.array(fixed_x)) 
 20              y_to_x[y] = scipy.array(fixed_x) 
 21          else: 
 22              for this_x in range_to_list(ranges[0]): 
 23                  fix_part(f,ranges[1:],fixed_x+[this_x]) 
 24      fix_part(f,ranges,[]) 
 25      allys = y_to_x.keys() 
 26      allys.sort() 
 27      return y_to_x[allys[0]] 
 28                   
 29   
 30 -def fmin_gen(f,x0,population_size=100,survivors=20,generations=20000,intruders=0,use_pp=True,convenience_class=None): 
 31      import struct 
 32      pmin = 1e-4 
 33      pmax = 1e4 
 34      def local_optimize(indiv,convenience_class=None): 
 35          better_indiv = indiv 
 36          if convenience_class: 
 37              better_indiv = scipy.optimize.fmin_bfgs(convenience_class.f,better_indiv,disp=0,maxiter=20) 
 38              fval = convenience_class.f(better_indiv) 
 39          else: 
 40              #better_indiv = scipy.optimize.fmin_bfgs(f,better_indiv,disp=0,maxiter=20) 
 41              fval = f(better_indiv) 
 42          return [better_indiv,fval] 
 43      def floatToBits(value): 
 44          if type(value) == type(1.0): 
 45              return (str(struct.unpack('Q', struct.pack('d', value))[0])).rjust(20,"0") 
 46          else: 
 47              return ",".join([floatToBits(x) for x in value.tolist()]) 
 48      def bitsToFloat(bits): 
 49          if "," not in bits: 
 50              return struct.unpack('d', struct.pack('Q', long(bits)))[0] 
 51          else: 
 52              return scipy.array([bitsToFloat(x) for x in bits.split(",")]) 
 53      def new_individual(): 
 54          #x = -scipy.log(scipy.array([scipy.rand() for i in range(indiv_size)])) 
 55          x = scipy.exp(scipy.rand(indiv_size)*10-5) 
 56          return x 
 57      def bool_mate(mother,father): 
 58          ms = floatToBits(mother) 
 59          fs = floatToBits(father) 
 60          l = random.sample(range(len(ms)),3) 
 61          l.sort() 
 62          [i1,i2,i3] = l 
 63          cs = ms[:i1]+fs[i1:i2]+ms[i2:i3]+fs[i3:] 
 64          child = bitsToFloat(cs) 
 65          if child.size != mother.size or None in child or scipy.inf in child or -scipy.inf in child: raise ValueError() 
 66          return child 
 67      def mate(mflist): 
 68          #return bool_mate(mflist[0],mflist[1]) 
 69          return mflist[0]+mflist[1]-mflist[2] 
 70      def limit(indiv): 
 71          for i in range(len(indiv.tolist())): 
 72              if indiv[i] > pmax: 
 73                  indiv[i] = pmax 
 74              elif indiv[i] < pmin: 
 75                  indiv[i] = pmin 
 76          return indiv 
 77      def mutate(indiv): 
 78          if not False: 
 79              bi = floatToBits(indiv) 
 80              number = int(math.ceil(float(len(bi))/100.)) 
 81              change_indices = random.sample([x for x in range(3,len(bi))],number) 
 82              for ci in change_indices: 
 83                  new = bi[:ci]+str(random.choice(range(10)))+bi[(ci+1):] 
 84                  try: 
 85                      bitsToFloat(new) 
 86                      bi = new 
 87                  except struct.error: 
 88                      # dont accept change 
 89                      pass 
 90              return scipy.absolute(bitsToFloat(bi)) 
 91          else: 
 92              return scipy.exp(scipy.log(indiv)+wolf["mutation_factor"]*scipy.array([random.normalvariate(0,1) for x in range(indiv_size)])) 
 93   
 94      if use_pp: 
 95          import pp 
 96          # get pp servers 
 97          servers = [] 
 98          fi = file("servers","r") 
 99          for l in fi: 
100              if l.startswith("#"): continue 
101              servers.append(l.split(",")[1]) 
102          fi.close() 
103          # read secret file 
104          try: 
105              fi = open("secret","r") 
106              secret = fi.readlines()[0] 
107              fi.close() 
108          except: 
109              import sys 
110              print "Please create a file calles \"secret\" with one long string in the first line." 
111              sys.exit(1) 
112   
113          print "starting servers",servers 
114          job_server = pp.Server(ppservers=tuple(servers),secret=secret) 
115   
116   
117   
118      population = [x0] 
119      indiv_size = x0.size 
120      quality = [] 
121      for i in range(population_size-len(population)): 
122          population.append(new_individual()) 
123      best_indiv = copy.deepcopy(x0) 
124      try: 
125          for i in range(generations): 
126              # rate individuals 
127              local_optimization_results = [] 
128              pre_computed_qualities = len(quality) 
129              for j in range(pre_computed_qualities,population_size): 
130                  if use_pp: 
131                      if convenience_class: 
132                          local_optimization_results.append( job_server.submit(local_optimize,(population[j],convenience_class),modules=("scipy","scipy.optimize","scipy.linalg","nnca"),globals=globals() ) ) 
133                      else: 
134                          local_optimization_results.append( job_server.submit(local_optimize,(population[j],),modules=("scipy","scipy.optimize","scipy.linalg","nnca"),globals=globals() ) ) 
135                  else: 
136                      local_optimization_results.append( local_optimize(population[j]) ) 
137              if use_pp: 
138                  for j in range(len(local_optimization_results)): 
139                      local_optimization_results[j] = local_optimization_results[j]() 
140              for j in range(len(local_optimization_results)): 
141                  [better_indiv,fval] = local_optimization_results[j] 
142                  population[pre_computed_qualities+j] = better_indiv 
143                  quality.append( fval ) 
144              # replace None results 
145              for j in range(len(population)): 
146                  while quality[j] == None: 
147                      population[j] = new_individual() 
148                      quality[j] = f(population[j]) 
149              # sort 
150              sorted_quality = zip(quality,population) 
151              sorted_quality.sort(key=lambda x:x[0]) 
152              new_population = [] 
153              new_quality = [] 
154              for j in range(survivors): 
155                  f_val,indiv = sorted_quality[j] 
156                  new_population.append(indiv) 
157                  new_quality.append(f_val) 
158              population = new_population 
159              quality = new_quality 
160              # intrude 
161              for j in range(intruders): 
162                  population.append(new_individual()) 
163              # mate 
164              current_size = len(population) 
165              for j in range(current_size,population_size): 
166                  population.append(limit(mutate(mate(random.sample(population[:current_size],3))))) 
167              print "generation",str(i+1).rjust(7),"     f =",quality[0] 
168              best_indiv = copy.deepcopy(population[0]) 
169              print best_indiv 
170      except KeyboardInterrupt: 
171          print "Stopping computation" 
172      if use_pp: 
173          job_server.print_stats() 
174      print "generation goodbye      f =",quality[0] 
175      return best_indiv 
176   
177 -def fmin_differential_evolution(f,x0,population_size=100,generations=20000,bounds=[1e-4,1e4],crossover_factor=0.2): 
178      import struct 
179      def floatToBits(value): 
180          if type(value) == type(1.0): 
181              return (str(struct.unpack('Q', struct.pack('d', value))[0])).rjust(20,"0") 
182          else: 
183              return ",".join([floatToBits(x) for x in value.tolist()]) 
184      def bitsToFloat(bits): 
185          if "," not in bits: 
186              return struct.unpack('d', struct.pack('Q', long(bits)))[0] 
187          else: 
188              return scipy.array([bitsToFloat(x) for x in bits.split(",")]) 
189      def new_individual(): 
190          logmin = scipy.log(bounds[0]) 
191          logmax = scipy.log(bounds[1]) 
192          x = scipy.exp(scipy.rand(indiv_size)*(logmax-logmin)-logmin) 
193          return x 
194      def crossover(orig,crossing_vector): 
195          for i in range(len(orig.tolist())): 
196              if random.random()<crossover_factor: 
197                  orig[i] = crossing_vector[i] 
198      population = [x0] 
199      indiv_size = x0.size 
200      quality = [f(x0)] 
201      for i in range(population_size-len(population)): 
202          population.append(new_individual()) 
203          quality.append(f(population[-1])) 
204      best_indiv = copy.deepcopy(x0) 
205      try: 
206          for i in range(generations): 
207              for targetindex in range(len(population)): 
208                  [v1,v2,v3] = random.sample(population,3) 
209                  trial_vector = v1+v2-v3 
210                  target_vector = population[targetindex] 
211                  crossover(trial_vector,target_vector) 
212                  trial_quality = f(trial_vector) 
213                  if trial_quality < quality[targetindex]: 
214                      population[targetindex] = trial_vector 
215                      quality[targetindex] = trial_quality 
216                   
217              # replace None results 
218              for j in range(len(population)): 
219                  while quality[j] == None: 
220                      population[j] = new_individual() 
221                      quality[j] = f(population[j]) 
222   
223              best_index = quality.index(min(quality)) 
224              print "generation",str(i+1).rjust(7),"     f =",quality[best_index] 
225              print population[best_index] 
226              best_indiv = copy.deepcopy(population[best_index]) 
227      except KeyboardInterrupt: 
228          print "Stopping computation" 
229      return best_indiv 
230   
231 -def make_spaces_in_formula(formula): 
232      f = " "+copy.deepcopy(formula)+" " 
233      for p in "()+-*/": 
234          f = f.replace(p," "+p+" ") 
235      f = f.replace("  "," ") 
236      return f 
237   
238 -class tf_function: 
239 -    def __init__(self,constant_part,constant_parameters,role2function_per_tf,tf_parameters,tf_delimiter,total_number_tf): 
240          """In the constant_part tf_part is replaced by the tf_function which is (function_per_tf) tf_delimiter (function_per_tf) and in function_per_tf tf describes the tf.""" 
241          self.constant_part = make_spaces_in_formula(constant_part) 
242          self.constant_parameters = constant_parameters 
243          self.role2function_per_tf = role2function_per_tf 
244          self.tf_parameters = tf_parameters 
245          self.tf_delimiter = tf_delimiter 
246          self.total_number_tf = total_number_tf 
247 -    def generate_functions_for_role(self,tf2roles): 
248          # generating the constant part 
249          self.all_parameters = copy.deepcopy(self.constant_parameters) 
250          self.complete_function = copy.deepcopy(self.constant_part) 
251          # generating the tf part 
252          functions_for_all_tfs = [] 
253          new_tf_names = [] 
254          for index in range(self.total_number_tf): 
255              # copy the parameters 
256              this_tf_function = make_spaces_in_formula(copy.deepcopy(self.role2function_per_tf[tf2roles[index]])) 
257              newtf = "tf_"+str(index) 
258              this_tf_function = this_tf_function.replace(" tf "," "+newtf+" ") 
259              new_tf_names.append(newtf) 
260              for p in self.tf_parameters: 
261                  newp = p+"_"+str(index) 
262                  this_tf_function = this_tf_function.replace(p,newp) 
263              self.all_parameters.append(newp) 
264              functions_for_all_tfs.append(this_tf_function) 
265          for ftn in new_tf_names: 
266              self.all_parameters.append(ftn) 
267          function_for_all_tfs = self.tf_delimiter.join(functions_for_all_tfs) 
268          self.complete_function = self.complete_function.replace(" x ",function_for_all_tfs) 
269 -    def __evaluate_function__(self,all_parameter_values): 
270          all_parameter_dict = dict(zip(self.all_parameters,all_parameter_values)) 
271          return eval(self.complete_function,all_parameter_dict) 
272 -    def evaluate_function(self,constant_parameter_values,tf_parameter_values,tf_values,tf2roles): 
273          self.generate_functions_for_role(tf2roles) 
274          return self.__evaluate_function__(constant_parameter_values.tolist()+tf_parameter_values.tolist()+tf_values.tolist()) 
275   
276 -class tf_system: 
277 -    def __init__(self,constant_part,constant_parameters,role2function_per_tf,tf_parameters,tf_delimiter,number_tf,protein_matrix,force_positive,role_matrix,role2free_param_index): 
278          self.zero_conditions_in_forcing_function = True 
279          self.protein_matrix = protein_matrix 
280          self.number_time_points = self.protein_matrix.shape[1] 
281          self.number_proteins = self.protein_matrix.shape[0] 
282          self.number_tf = number_tf 
283          self.force_positive = force_positive 
284          self.tf_function = tf_function(constant_part,constant_parameters,role2function_per_tf,tf_parameters,tf_delimiter,self.number_tf) 
285          self.number_constant_parameters = len(constant_parameters) 
286          self.number_tf_parameters = len(tf_parameters) 
287          self.number_total_parameters = self.number_proteins * len(constant_parameters) + self.number_proteins * self.number_tf * len(tf_parameters) + self.number_time_points * self.number_tf 
288          self.x = scipy.array([1.0]*self.number_total_parameters) 
289          self.role_matrix = role_matrix 
290          self.role2free_param_index = role2free_param_index 
291          self.force_zero_indices = self.get_force_zero_indices() 
292 -    def evaluate_function(self,x): 
293          group2offset = self.number_constant_parameters * self.number_proteins 
294          group3offset = group2offset + self.number_proteins * self.number_tf * self.number_tf_parameters 
295          sum_of_squares = 0.0 
296          for protein_index in range(self.number_proteins): 
297              for time_index in range(self.number_time_points): 
298                  target_value = self.protein_matrix[protein_index,time_index] 
299                  # getting the constant parameters 
300                  start_index = protein_index*self.number_constant_parameters 
301                  stop_index = start_index+self.number_constant_parameters 
302                  constant_parameter_values = x[start_index:stop_index] 
303                  # getting the constant parameters 
304                  start_index = group2offset + protein_index * self.number_tf * self.number_tf_parameters 
305                  stop_index = start_index+self.number_tf * self.number_tf_parameters 
306                  tf_parameter_values = x[start_index:stop_index] 
307                  # getting the constant parameters 
308                  start_index = group3offset+time_index*self.number_tf 
309                  stop_index = start_index+self.number_tf 
310                  tf_values = x[start_index:stop_index] 
311                  # evaluate 
312                  estimated_value = self.tf_function.evaluate_function(constant_parameter_values,tf_parameter_values,tf_values,self.role_matrix[protein_index]) 
313                  sum_of_squares += (target_value-estimated_value)**2 
314          if self.force_positive: 
315              sum_of_squares += scipy.linalg.norm(x-scipy.absolute(x)) 
316          if self.zero_conditions_in_forcing_function: 
317              for index in self.force_zero_indices: 
318                  sum_of_squares += x[index]**2 
319          # scaling everything 
320          # sum_of_squares += (max(x[group3offset:])-1)**2 
321          return sum_of_squares 
322 -    def f(self,x): 
323          return self.evaluate_function(x) 
324 -    def get_force_zero_indices(self): 
325          # penalty for free parameters which are not 0 
326          indices = [] 
327          for protein_index in range(self.number_proteins): 
328              for tf_index in range(self.number_tf): 
329                  thisrole = self.role_matrix[protein_index][tf_index] 
330                  free_param_indices = self.role2free_param_index[thisrole] 
331                  for free_param_index in free_param_indices: 
332                      index_in_x = self.number_constant_parameters * self.number_proteins + protein_index * self.number_tf * self.number_tf_parameters + tf_index * self.number_tf_parameters + free_param_index 
333                      indices.append(index_in_x) 
334          return indices 
335 -    def __repr__(self): 
336          s  = "Transcription factor system\n" 
337          s += "--- Data: v Protein > Time\n" 
338          s += str(self.protein_matrix)+"\n" 
339          if self.number_constant_parameters>0: 
340              s += "--- Constant parameters: v Protein > Parameter\n" 
341              m = scipy.matrix([[self.x[pi*self.number_constant_parameters+ci] for ci in range(self.number_constant_parameters)] for pi in range(self.number_proteins)]) 
342              s += str(m)+"\n" 
343          if self.number_tf_parameters>0: 
344              s += "--- TF  parameters: v Protein > TF [> Parameter]\n" 
345              for tf_param_index in range(self.number_tf_parameters): 
346                  s += "- Param "+str(self.tf_function.tf_parameters[tf_param_index])+"\n" 
347                  m = scipy.matrix([[self.x[self.number_constant_parameters*self.number_proteins+pi*self.number_tf*self.number_tf_parameters+tfi*self.number_tf_parameters+tf_param_index] for tfi in range(self.number_tf)] for pi in range(self.number_proteins)]) 
348                  s += str(m)+"\n" 
349          s += "--- TF concentrations: v TF > Time\n" 
350          g3o = self.number_constant_parameters * self.number_proteins + self.number_proteins * self.number_tf * self.number_tf_parameters 
351          m = scipy.matrix([[self.x[g3o+ti*self.number_tf+tfi] for ti in range(self.number_time_points)] for tfi in range(self.number_tf)]) 
352          s += str(m)+"\n" 
353          s += "Objective value: "+str(self.evaluate_function(self.x))+"\n" 
354          return s 
355 -    def __str__(self): 
356          return repr(self) 
357 -    def estimate(self): 
358          def f(x): 
359              return self.evaluate_function(x) 
360          self.x = fmin_gen(f,self.x,convenience_class=self) 
361          #self.x = scipy.optimize.fmin_bfgs(f,self.x) 
362          for index in self.force_zero_indices: 
363              self.x[index] = 0. 
364          save_tfs(self) 
365   
366                   
367           
368 -class dataloader: 
369 -    def __init__(self,filename,delimiter): 
370          self.lines = [] 
371          f = open(filename,"r") 
372          for line in f: 
373              self.lines.append([float(x) for x in line.strip("\n").split(delimiter)]) 
374          f.close() 
375 -    def get_matrix(self): 
376          return scipy.matrix(self.lines) 
377   
378 -def load_tfs(): 
379      try: 
380          f = file("tfsystem","r") 
381      except IOError: 
382          return None 
383      tfs = pickle.load(f) 
384      f.close() 
385      return tfs 
386 -def save_tfs(tfs): 
387      f = file("tfsystem","w") 
388      pickle.dump(tfs,f) 
389      f.close() 
390   
391 -def main(): 
392      constant_part = "c*x" 
393      constant_parameters = ["c"] 
394      role2function_per_tf = {1:"(tf/a)/(1+tf/a)",0:"0"} 
395      tf_parameters = ["a"] 
396      tf_delimiter = "*" 
397      number_of_tf = 2 
398      #protein_matrix = scipy.matrix([[1.,2.,3.,4.,5.],[2.,3.,4.,5.,6.],[3.,5.,7.,9.,11.],[6.,10.,14.,18.,22.]]) 
399      protein_matrix = dataloader("spassdaten",";").get_matrix() 
400      role_matrix = [[1,1],[1,1],[1,1],[1,1]] 
401      role2free_param_index = {1:[],0:[0]} 
402      tfs = load_tfs() 
403      if not tfs: 
404          tfs = tf_system(constant_part,constant_parameters,role2function_per_tf,tf_parameters,tf_delimiter,number_of_tf,protein_matrix,True,role_matrix,role2free_param_index) 
405      tfs.estimate() 
406      print tfs 
407   
408   
409   
410  if __name__ == "__main__": 
411      #main() 
412 -    def f(x): 
413          for i in range(len(x.tolist())): 
414              if abs(x[i])>10e3: 
415                  return abs(x[i]) 
416          return (x[0] - 4)**2 + (x[1]-8)**2 + (x[2]-1)**2 + (x[3]-9)**2 + (x[4]-8)**2 + (x[5]-2)**2 
417       
418      fmin_differential_evolution(f,scipy.array([1,1,1,1,1,1]),10000,20000,[1e-1,1e1],0.2) 
419