Source Code for Module model_selection.bootstrap

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