Source Code for Module model_selection.exp_overall_disc

  1  # Model selection using the expected overall discrepancy. 
  2  # 
  3  # The input relaxation data for this method should be the true data (theoretical, back calculated relaxation values). 
  4  # 
  5  # The program is divided into the following stages: 
  6  #       Stage 1:  Creation of the files for the model-free calculations for models 1 to 5.  Monte Carlo 
  7  #               simulations are used, but the initial data rather than the backcalculated data is randomised. 
  8  #       Stage 2:  Model selection and the creation of the final run.  Monte Carlo simulations are used to 
  9  #               find errors.  This stage has the option of optimizing the diffusion tensor along with the 
 10  #               model-free parameters. 
 11  #       Stage 3:  Extraction of the data. 
 12   
 13  import sys 
 14  from math import log, pi 
 15  from re import match 
 16   
 17  from common_ops import common_operations 
 18   
 19   
 20 -class exp_overall_disc(common_operations): 
 21 -        def __init__(self, mf): 
 22                  "Model-free analysis based on the expected overall discrepancy." 
 23   
 24                  self.mf = mf 
 25   
 26   
 27 -        def model_selection(self): 
 28                  "Model selection." 
 29   
 30                  data = self.mf.data.data 
 31                  self.mf.data.calc_frq() 
 32                  self.mf.data.calc_constants() 
 33                  n = float(self.mf.data.num_data_sets) 
 34                  tm = float(self.mf.usr_param.tm['val']) * 1e-9 
 35   
 36                  if self.mf.debug: 
 37                          self.mf.log.write("\n\n<<< Expected overall discrepancy >>>\n\n") 
 38   
 39                  print "Calculating the expected overall discrepancy" 
 40                  for res in range(len(self.mf.data.relax_data[0])): 
 41                          print "Residue: " + self.mf.data.relax_data[0][res][1] + " " + self.mf.data.relax_data[0][res][0] 
 42                          self.mf.data.results.append({}) 
 43                          file_name = self.mf.data.relax_data[0][res][1] + '_' + self.mf.data.relax_data[0][res][0] + '.out' 
 44   
 45                          if self.mf.debug: 
 46                                  self.mf.log.write('%-22s\n' % ( "< Checking res " + data['m1'][res]['res_num'] + " >\n")) 
 47   
 48                          real = [] 
 49                          err = [] 
 50                          types = [] 
 51                          for i in range(self.mf.data.num_ri): 
 52                                  real.append(float(self.mf.data.relax_data[i][res][2])) 
 53                                  err.append(float(self.mf.data.relax_data[i][res][3])) 
 54                                  types.append([self.mf.data.data_types[i], float(self.mf.data.frq[self.mf.data.remap_table[i]])]) 
 55   
 56                          for model in self.mf.data.runs: 
 57                                  if self.mf.debug: 
 58                                          self.mf.log.write("\nCalculating expected overall discrepancy for res " + `res` + ", model " + model + "\n\n") 
 59                                          for i in range(self.mf.data.num_ri): 
 60                                                  self.mf.log.write("-------------------") 
 61                                          self.mf.log.write("\n") 
 62                                          for i in range(self.mf.data.num_ri): 
 63                                                  name = " Orig " + self.mf.data.frq_label[self.mf.data.remap_table[i]] + " " + self.mf.data.data_types[i] 
 64                                                  self.mf.log.write("%-17s%2s" % (name, " |")) 
 65                                          self.mf.log.write("\n") 
 66                                          for i in range(self.mf.data.num_ri): 
 67                                                  self.mf.log.write("%8.4f" % self.mf.data.relax_data[i][res][2]) 
 68                                                  self.mf.log.write("%1s" % "±") 
 69                                                  self.mf.log.write("%-8.4f" % self.mf.data.relax_data[i][res][3]) 
 70                                                  self.mf.log.write("%2s" % " |") 
 71                                          self.mf.log.write("\n") 
 72                                          for i in range(self.mf.data.num_ri): 
 73                                                  self.mf.log.write("-------------------") 
 74                                          self.mf.log.write("\n\n") 
 75   
 76                                  file = self.mf.file_ops.open_file(model + "/" + file_name) 
 77                                  sum_chi2 = 0.0 
 78                                  num_sims = float(len(file)) 
 79                                  for sim in range(len(file)): 
 80                                          if self.mf.debug: 
 81                                                  self.mf.log.write("%5s%-10i%2s" % ("Sim: ", sim, " |")) 
 82   
 83                                          if match('m1', model): 
 84                                                  back_calc = self.mf.calc_relax_data.calc(tm, model, types, [ file[sim][2] ]) 
 85                                          elif match('m2', model) or match('m3', model): 
 86                                                  back_calc = self.mf.calc_relax_data.calc(tm, model, types, [ file[sim][2], file[sim][3] ]) 
 87                                          elif match('m4', model) or match('m5', model): 
 88                                                  back_calc = self.mf.calc_relax_data.calc(tm, model, types, [ file[sim][2], file[sim][3], file[sim][4] ]) 
 89                                          chi2 = self.mf.calc_chi2.relax_data(real, err, back_calc) 
 90                                          sum_chi2 = sum_chi2 + chi2 
 91   
 92                                          if self.mf.debug: 
 93                                                  self.mf.log.write("%7s%-10.4f%2s" % (" Chi2: ", chi2, " |")) 
 94                                                  self.mf.log.write("%11s%-13.4f%2s\n" % (" Sum Chi2: ", sum_chi2, " |")) 
 95   
 96                                  ave_chi2 = sum_chi2 / num_sims 
 97   
 98                                  if self.mf.debug: 
 99                                          self.mf.log.write("\nAverage Chi2 is: " + `ave_chi2` + "\n\n") 
100   
101                                  data[model][res]['expect'] = ave_chi2 / (2.0 * n) 
102   
103                          # Select model. 
104                          min = 'm1' 
105                          for model in self.mf.data.runs: 
106                                  if data[model][res]['expect'] < data[min][res]['expect']: 
107                                          min = model 
108                          if data[min][res]['expect'] == float('inf'): 
109                                  self.mf.data.results[res] = self.fill_results(data[min][res], model='0') 
110                          else: 
111                                  self.mf.data.results[res] = self.fill_results(data[min][res], model=min[1]) 
112   
113                          if self.mf.debug: 
114                                  self.mf.log.write(self.mf.usr_param.method + " (m1): " + `data['m1'][res]['expect']` + "\n") 
115                                  self.mf.log.write(self.mf.usr_param.method + " (m2): " + `data['m2'][res]['expect']` + "\n") 
116                                  self.mf.log.write(self.mf.usr_param.method + " (m3): " + `data['m3'][res]['expect']` + "\n") 
117                                  self.mf.log.write(self.mf.usr_param.method + " (m4): " + `data['m4'][res]['expect']` + "\n") 
118                                  self.mf.log.write(self.mf.usr_param.method + " (m5): " + `data['m5'][res]['expect']` + "\n") 
119                                  self.mf.log.write("The selected model is: " + min + "\n\n") 
120   
121                          print "   Model " + self.mf.data.results[res]['model'] 
122   
123   
124 -        def print_data(self): 
125                  "Print all the data into the 'data_all' file." 
126   
127                  file = open('data_all', 'w') 
128                  file_crit = open('crit', 'w') 
129   
130                  sys.stdout.write("[") 
131                  for res in range(len(self.mf.data.results)): 
132                          sys.stdout.write("-") 
133                          file.write("\n\n<<< Residue " + self.mf.data.results[res]['res_num']) 
134                          file.write(", Model " + self.mf.data.results[res]['model'] + " >>>\n") 
135                          file.write('%-20s' % '') 
136                          file.write('%-19s' % 'Model 1') 
137                          file.write('%-19s' % 'Model 2') 
138                          file.write('%-19s' % 'Model 3') 
139                          file.write('%-19s' % 'Model 4') 
140                          file.write('%-19s' % 'Model 5') 
141   
142                          file_crit.write('%-6s' % self.mf.data.results[res]['res_num']) 
143                          file_crit.write('%-6s' % self.mf.data.results[res]['model']) 
144   
145                          # S2. 
146                          file.write('\n%-20s' % 'S2') 
147                          for model in self.mf.data.runs: 
148                                  file.write('%9.3f' % self.mf.data.data[model][res]['s2']) 
149                                  file.write('%1s' % '±') 
150                                  file.write('%-9.3f' % self.mf.data.data[model][res]['s2_err']) 
151   
152                          # S2f. 
153                          file.write('\n%-20s' % 'S2f') 
154                          for model in self.mf.data.runs: 
155                                  file.write('%9.3f' % self.mf.data.data[model][res]['s2f']) 
156                                  file.write('%1s' % '±') 
157                                  file.write('%-9.3f' % self.mf.data.data[model][res]['s2f_err']) 
158   
159                          # S2s. 
160                          file.write('\n%-20s' % 'S2s') 
161                          for model in self.mf.data.runs: 
162                                  file.write('%9.3f' % self.mf.data.data[model][res]['s2s']) 
163                                  file.write('%1s' % '±') 
164                                  file.write('%-9.3f' % self.mf.data.data[model][res]['s2s_err']) 
165   
166                          # te. 
167                          file.write('\n%-20s' % 'te') 
168                          for model in self.mf.data.runs: 
169                                  file.write('%9.2f' % self.mf.data.data[model][res]['te']) 
170                                  file.write('%1s' % '±') 
171                                  file.write('%-9.2f' % self.mf.data.data[model][res]['te_err']) 
172   
173                          # Rex. 
174                          file.write('\n%-20s' % 'Rex') 
175                          for model in self.mf.data.runs: 
176                                  file.write('%9.3f' % self.mf.data.data[model][res]['rex']) 
177                                  file.write('%1s' % '±') 
178                                  file.write('%-9.3f' % self.mf.data.data[model][res]['rex_err']) 
179   
180                          # Chi2. 
181                          file.write('\n%-20s' % 'Chi2') 
182                          for model in self.mf.data.runs: 
183                                  file.write('%-19.3f' % self.mf.data.data[model][res]['chi2']) 
184   
185                          # Expected overall discrepancy. 
186                          file.write('\n%-20s' % 'Expect') 
187                          for model in self.mf.data.runs: 
188                                  file.write('%-19.3f' % self.mf.data.data[model][res]['expect']) 
189   
190                                  file_crit.write('%-25s' % `self.mf.data.data[model][res]['expect']`) 
191                          file_crit.write('\n') 
192   
193                  file.write('\n') 
194                  sys.stdout.write("]\n") 
195                  file.close() 
196