Source Code for Module modelmage.discriminator

  1  #-----------------------------------------------# 
  2  #                  modelMaGe                    # 
  3  #                                               # 
  4  #                 v1.0beta September 2009           # 
  5  #       Joerg Schaber, Max Floettmann and Jian Li   # 
  6  #                                               # 
  7  #                http://modelmage.org               # 
  8  #-----------------------------------------------# 
  9  ## @file    discriminator.py 
 10  ## @brief   Provides functions for comparison and output of copasi result files. 
 11  ## @author  Max Floettmann and Jian Li 
 12  ##  
 13  ## This file is part of modelMaGe.  Please visit http://modelMaGe.org for more 
 14  ## information about modelMaGe, and the latest version of modelMaGe. 
 15  ## 
 16   
 17  import os 
 18  import re 
 19  import glob 
 20  import math 
 21   
 22   
 23 -class Discriminator: 
 24          """ 
 25          Provides functions for comparison and output of copasi result files. 
 26          """ 
 27   
 28 -        def __init__(self, ResultCopasi): 
 29                  ''' 
 30                  @type ResultCopasi: str 
 31                  @param ResultCopasi: path to the created copasifile              
 32                   
 33                  to get the parameter estimation result file at the right place 
 34                  ''' 
 35                   
 36                  self.ResultCopasi = ResultCopasi 
 37   
 38 -        def _getObjectiveValue(self, result): 
 39                  """ 
 40                  Looks up the objective value in a Copasi result file and returns it. 
 41                  @type result: str 
 42                  @param result: path to the result file 
 43   
 44                  @rtype: str/bool 
 45                  @return: returns the value if it is found, otherwise returns False 
 46                  """ 
 47                  p = re.compile('Objective Function Value:\s*(\S*)') 
 48                  m = re.findall(p, result) 
 49                  if m: 
 50                          m = m[-1] 
 51                          return m 
 52                  else: 
 53                          return False 
 54           
 55 -        def _getLastResult(self, file): 
 56                  f = open(file) 
 57                  text = f.read() 
 58                  p = re.compile("Parameter Estimation Result:") 
 59                  resultlist = re.split(p, text) 
 60                  newestResult = resultlist[-1] 
 61                  return newestResult 
 62           
 63 -        def _getNumberOfParameters(self, result): 
 64                  p = re.findall("(?:.*)\..*:\s*.*", result) # matches the parameter lines 
 65                  return len(p) 
 66           
 67 -        def _getNumberOfObservations(self, result): 
 68  #               experimentName = re.findall("Experiment:\s(.*)", result) 
 69                  data = re.findall("(?m)Row\s.*?$[^ A-Z]*", result) 
 70                  dataPoints = 0 
 71                  for experiment in data: 
 72                          obs = re.findall("(\[[\S]*\]\(Data\))", experiment) 
 73                          lines = re.findall("\d{1,5}\.\t", experiment) 
 74                          dataPoints += len(lines)*len(obs) 
 75           
 76                  return dataPoints 
 77   
 78 -        def AIC(self, rss, k, n): 
 79                  """ 
 80                  Calculates the Akaike information coefficient to get a better ranking for 
 81                  the models. 
 82                  """ 
 83                  aic = 2*k + n*((math.log(2*math.pi*float(rss)/n))+1) 
 84                  return aic 
 85                   
 86 -        def AICc(self, aic, n, k): 
 87                  if(n>k-1): 
 88                          aicc = aic + 2*((k*(k+1))/(n-k-1)) 
 89                  else: 
 90                          aicc = 'nan' 
 91                  return aicc 
 92   
 93 -        def getRanking(self): 
 94                  """ 
 95                  Prints a table of the files in './ResultCopasiFiles' and their objective values. 
 96                   
 97                  '#\tModel\tObjective Value' 
 98   
 99                  @rtype: None 
100                  @return: none 
101                  """ 
102                  ranking = [] 
103                  longestName = 0 
104                  # __init__() method is only for this purpose 
105                  tmp = os.path.join(os.path.abspath(os.path.curdir), self.ResultCopasi) 
106                  if not os.path.exists(tmp): 
107                          print "Could not find the Resultfiles in %s." % tmp 
108                  for file in glob.glob1(tmp,"*_est.txt"): 
109                          file = os.path.join(tmp,file) 
110                          lastResult = self._getLastResult(file) 
111                          rss = self._getObjectiveValue(lastResult) 
112                          n = self._getNumberOfObservations(lastResult) 
113                          k = self._getNumberOfParameters(lastResult) 
114  #                       print "file: %s\t rss: %s\tn: %s\tk: %s"%(file,rss, n, k) 
115                          try: 
116                                  aic = self.AIC(rss, k, n) 
117                                  aicc = self.AICc(aic, n, k) 
118                          except OverflowError: 
119                                  print 'No standard report found in reportfile. Please reset report settings in Copasi and start again.' 
120                                  exit(1) 
121                          except ZeroDivisionError: 
122                                  print "Could not find standard report elements in the result file %s. Exiting ..." % file 
123                                  exit(1) 
124                          ranking.append((file,n,k,rss,aic,aicc)) 
125                          filename = os.path.split(file)[1] 
126                          if len(filename) > longestName: 
127                                  longestName = len(filename) 
128                  ranking.sort(key=lambda x: x[5]) 
129                  rank = 1 
130                  out = '' 
131                  out += '\nRanking of models by AICc:\n' 
132                  out += '%3s %s%5s%5s%20s%10s'% ('#', 'Model'.rjust(longestName+5),'n','k','Objective Value','AICc\n') 
133                  for elem in ranking: 
134                          filename = os.path.split(elem[0])[1] 
135                          out += '%3d.%s%5d%5d%20.4f%10.4f\n' % (rank, filename.rjust(longestName+5), elem[1], elem[2], float(elem[3]), float(elem[5])) 
136                          rank += 1 
137                  return out 
138   
139   
140           
141 -def main(): 
142          myPath = '/result/ResultCopasiFiles'  
143          d = Discriminator(myPath) 
144          print d.getRanking() 
145   
146   
147  if __name__ == '__main__':       
148          main() 
149