Source Code for Module true

  1  # Calculate the overall discrepency. 
  2  # 
  3  # The input relaxation data for this method should be the operating model data (the theoretical, back calculated 
  4  # relaxation values) which has been Gaussian randomized for a given error.  The original operating model data 
  5  # should be placed in the file 'op_data'.  The format of the file should be as follows.  The first line is a 
  6  # header line and is ignored.  The columns are: 
  7  #       0 - Residue number. 
  8  #       1 - Residue name. 
  9  #       2 - R1 operating model data (1st field strength). 
 10  #       3 - R2 operating model data (1st field strength). 
 11  #       4 - NOE operating model data (1st field strength). 
 12  #       5 - R1 operating model data (2nd field strength). 
 13  #       6 - R2 operating model data (2nd field strength). 
 14  #       7 - NOE operating model data (2nd field strength). 
 15  #       8 - etc 
 16  # 
 17  # 
 18  # The program is divided into the following stages: 
 19  #       Stage 1:   Creation of the files for the modelfree calculations for models 1 to 5 for the randomized 
 20  #               data. 
 21  #       Stage 2a:  Calculation of the overall discrepency for model selection, and the creation of the final 
 22  #               run.  Monte Carlo simulations are used to find errors, and the diffusion tensor is unoptimized. 
 23  #               Files are placed in the directory 'final'. 
 24  #       Stage 2b:  Calculation of the overall discrepency for model selection, and the creation of the final 
 25  #               optimization run.  Monte Carlo simulations are used to find errors, and the diffusion tensor 
 26  #               is optimized.  Files are placed in the directory 'optimized'. 
 27  #       Stage 3:   Extraction of optimized data. 
 28   
 29  import sys 
 30  from math import log, pi 
 31  from re import match 
 32   
 33  from common_ops import common_operations 
 34   
 35   
 36 -class true(common_operations): 
 37 -        def __init__(self, mf): 
 38                  "Calculation of the theoretical overall discrepency." 
 39   
 40                  self.mf = mf 
 41   
 42                  print "Modelfree analysis based on the overall discrepency for model selection." 
 43                  self.initialize() 
 44                  message = "See the file 'true.py' for details." 
 45                  self.mf.file_ops.read_file('op_data', message) 
 46                  self.mf.data.true.op_data = self.mf.file_ops.open_file(file_name='op_data') 
 47                  self.mf.data.runs = ['m1', 'm2', 'm3', 'm4', 'm5'] 
 48                  self.goto_stage() 
 49   
 50   
 51 -        def calc_crit(self, res, n, k, chisq): 
 52                  "Calculate the criteria" 
 53   
 54                  sum_ln_err = 0 
 55                  for i in range(len(self.mf.data.relax_data)): 
 56                          ln_err = log(float(self.mf.data.relax_data[i][res][3])) 
 57                          sum_ln_err = sum_ln_err + ln_err 
 58   
 59                  if match('^AIC$', self.mf.data.usr_param.method): 
 60                          aic = n*log(2*pi) + sum_ln_err + chisq + 2*k 
 61                          aic = aic / (2*n) 
 62                          return aic 
 63   
 64                  elif match('^AICc$', self.mf.data.usr_param.method): 
 65                          aicc = n*log(2*pi) + sum_ln_err + chisq + 2*k + 2*k*(k+1)/(n-k-1) 
 66                          aicc = aicc / (2*n) 
 67                          return aicc 
 68   
 69                  elif match('^BIC$', self.mf.data.usr_param.method): 
 70                          bic = n*log(2*pi) + sum_ln_err + chisq + k*log(n) 
 71                          bic = bic / ( 2 * n ) 
 72                          return bic 
 73   
 74   
 75 -        def model_selection(self): 
 76                  "Model selection." 
 77                   
 78                  data = self.mf.data.data 
 79   
 80                  self.mf.log.write("\n\n<<< " + self.mf.data.usr_param.method + " model selection >>>") 
 81                  for res in range(len(self.mf.data.relax_data[0])): 
 82                          self.mf.data.results.append({}) 
 83                          self.mf.log.write('\n%-22s' % ( "   Checking res " + data['m1'][res]['res_num'] )) 
 84   
 85                          n = self.mf.data.num_data_sets 
 86   
 87                          data['m1'][res]['crit'] = self.calc_crit(res, n, k=1, chisq=data['m1'][res]['sse']) 
 88                          data['m2'][res]['crit'] = self.calc_crit(res, n, k=2, chisq=data['m2'][res]['sse']) 
 89                          data['m3'][res]['crit'] = self.calc_crit(res, n, k=2, chisq=data['m3'][res]['sse']) 
 90                          data['m4'][res]['crit'] = self.calc_crit(res, n, k=3, chisq=data['m4'][res]['sse']) 
 91                          data['m5'][res]['crit'] = self.calc_crit(res, n, k=3, chisq=data['m5'][res]['sse']) 
 92   
 93                          # Select model. 
 94                          min = 'm1' 
 95                          for run in self.mf.data.runs: 
 96                                  if data[run][res]['crit'] < data[min][res]['crit']: 
 97                                          min = run 
 98                          self.mf.data.results[res] = self.fill_results(data[min][res], model=min[1]) 
 99   
100                          self.mf.log.write("\n\t" + self.mf.data.usr_param.method + " (m1): " + `data['m1'][res]['crit']` + "\n") 
101                          self.mf.log.write("\n\t" + self.mf.data.usr_param.method + " (m2): " + `data['m2'][res]['crit']` + "\n") 
102                          self.mf.log.write("\n\t" + self.mf.data.usr_param.method + " (m3): " + `data['m3'][res]['crit']` + "\n") 
103                          self.mf.log.write("\n\t" + self.mf.data.usr_param.method + " (m4): " + `data['m4'][res]['crit']` + "\n") 
104                          self.mf.log.write("\n\t" + self.mf.data.usr_param.method + " (m5): " + `data['m5'][res]['crit']` + "\n") 
105                          self.mf.log.write("\tThe selected model is: " + min + "\n\n") 
106   
107   
108 -        def stage2(self): 
109                  self.mf.file_ops.mkdir('grace') 
110                   
111                  print "\n[ Modelfree data extraction ]\n" 
112                  for run in self.mf.data.runs: 
113                          mfout = self.mf.file_ops.read_file(run + '/mfout') 
114                          mfout_lines = mfout.readlines() 
115                          mfout.close() 
116                          print "Extracting modelfree data from " + run + "/mfout." 
117                          num_res = len(self.mf.data.relax_data[0]) 
118                          self.mf.data.data[run] = self.mf.star.extract(mfout_lines, num_res) 
119   
120                  print "\n[ " + self.mf.data.usr_param.method + " model selection ]\n" 
121                  self.model_selection() 
122   
123                  print "\n[ Printing results ]\n" 
124                  self.print_results() 
125   
126                  print "\n[ Placing data structures into \"data_all\" ]\n" 
127                  self.print_data() 
128   
129                  print "\n[ Grace file creation ]\n" 
130                  self.grace('grace/S2.agr', 'S2', subtitle="After model selection, unoptimized") 
131                  self.grace('grace/S2f.agr', 'S2f', subtitle="After model selection, unoptimized") 
132                  self.grace('grace/S2s.agr', 'S2s', subtitle="After model selection, unoptimized") 
133                  self.grace('grace/te.agr', 'te', subtitle="After model selection, unoptimized") 
134                  self.grace('grace/Rex.agr', 'Rex', subtitle="After model selection, unoptimized") 
135                  self.grace('grace/SSE.agr', 'SSE', subtitle="After model selection, unoptimized") 
136