Source Code for Module model_selection.palmer

  1  # The method given by (Mandel et al., 1995) is divided into the three stages: 
  2  # 
  3  #       Stage 1:  Creation of the files for the initial model-free calculations for models 1 to 5, 
  4  #               and f-tests between them. 
  5  #       Stage 2:  Model selection and the creation of the final run.  Monte Carlo simulations are used to 
  6  #               find errors.  This stage has the option of optimizing the diffusion tensor along with the 
  7  #               model-free parameters. 
  8  #       Stage 3:  Extraction of the data. 
  9  # 
 10  # These stages are repeated until the data converges. 
 11   
 12  import sys 
 13  from re import match 
 14   
 15  from common_ops import common_operations 
 16   
 17   
 18 -class palmer(common_operations): 
 19 -        def __init__(self, mf): 
 20                  """The model-free analysis of Palmer. 
 21   
 22                  print "Palmer's method for model-free analysis. (Mandel et al., 1995)" 
 23                  """ 
 24                  self.mf = mf 
 25   
 26   
 27 -        def anova_tests(self): 
 28                  "Do the chi squared and F-tests." 
 29   
 30                  for res in range(len(self.mf.data.relax_data[0])): 
 31                          for run in self.mf.data.runs: 
 32                                  if match('^m', run): 
 33                                          # Chi squared test. 
 34                                          if self.mf.data.data[run][res]['chi2'] <= self.mf.data.data[run][res]['chi2_lim']: 
 35                                                  self.mf.data.data[run][res]['chi2_test'] = 1 
 36                                          else: 
 37                                                  self.mf.data.data[run][res]['chi2_test'] = 0 
 38   
 39                                          # Large chi squared test. 
 40                                          if self.mf.data.data[run][res]['chi2'] >= float(self.mf.usr_param.large_chi2): 
 41                                                  self.mf.data.data[run][res]['large_chi2'] = 1 
 42                                          else: 
 43                                                  self.mf.data.data[run][res]['large_chi2'] = 0 
 44   
 45                                          # Zero chi squared test. 
 46                                          if self.mf.data.data[run][res]['chi2'] == 0.0: 
 47                                                  self.mf.data.data[run][res]['zero_chi2'] = 1 
 48                                          else: 
 49                                                  self.mf.data.data[run][res]['zero_chi2'] = 0 
 50   
 51                                  else: 
 52                                          # F-test. 
 53                                          if self.mf.data.data[run][res]['fstat_lim'] < 1.5: 
 54                                                  if self.mf.data.data[run][res]['fstat'] > 1.5: 
 55                                                          self.mf.data.data[run][res]['ftest'] = 1 
 56                                                  else: 
 57                                                          self.mf.data.data[run][res]['ftest'] = 0 
 58                                          elif self.mf.data.data[run][res]['fstat_lim'] >= 1.5: 
 59                                                  if self.mf.data.data[run][res]['fstat'] > self.mf.data.data[run][res]['fstat_lim']: 
 60                                                          self.mf.data.data[run][res]['ftest'] = 1 
 61                                                  else: 
 62                                                          self.mf.data.data[run][res]['ftest'] = 0 
 63   
 64 -        def model_selection(self): 
 65                  "Palmer's model selection." 
 66   
 67                  self.anova_tests() 
 68                  if self.mf.data.num_data_sets > 3: 
 69                          # ie degrees of freedom > 0 in all models. 
 70                          print "The number of data sets is greater than 3." 
 71                          print "\tRunning Palmer's model selection, with additional chi-squared and F-tests" 
 72                          print "\tfor models 4 and 5 (the degrees of freedom for these models are greater than 0).\n" 
 73   
 74                          if self.mf.debug: 
 75                                  self.mf.log.write("Extended model selection.\n\n") 
 76   
 77                          self.model_selection_extended() 
 78                  else: 
 79                          print "The number of data sets is equal to 3." 
 80                          print "\tRunning Palmer's model selection.\n" 
 81   
 82                          if self.mf.debug: 
 83                                  self.mf.log.write("Normal model selection.\n\n") 
 84   
 85                          self.model_selection_normal() 
 86   
 87   
 88 -        def model_selection_normal(self): 
 89                  "Palmer's model selection." 
 90   
 91                  data = self.mf.data.data 
 92   
 93                  if self.mf.debug: 
 94                          self.mf.log.write("\n\n<<< Palmer's model selection >>>\n\n") 
 95   
 96                  for res in range(len(self.mf.data.relax_data[0])): 
 97                          self.mf.data.results.append({}) 
 98   
 99                          if self.mf.debug: 
100                                  self.mf.log.write('%-22s\n' % ( "Checking res " + data['m1'][res]['res_num'] )) 
101   
102                          # Model 1 test. 
103                          if data['m1'][res]['chi2_test']: 
104                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='1') 
105                                  if self.mf.debug: 
106                                          self.mf.log.write('%-12s' % '[Model 1]') 
107   
108                          # Test if both model 2 and 3 fit!!! (Should not occur) 
109                          elif data['m2'][res]['chi2_test'] and data['f-m1m2'][res]['ftest'] \ 
110                                  and data['m3'][res]['chi2_test'] and data['f-m1m3'][res]['ftest']: 
111                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='2+3') 
112                                  if self.mf.debug: 
113                                          self.mf.log.write('%-12s' % '[Model 2 and 3]') 
114   
115                          # Model 2 test. 
116                          elif data['m2'][res]['chi2_test'] and data['f-m1m2'][res]['ftest']: 
117                                  self.mf.data.results[res] = self.fill_results(data['m2'][res], model='2') 
118                                  if self.mf.debug: 
119                                          self.mf.log.write('%-12s' % '[Model 2]') 
120   
121                          # Model 3 test. 
122                          elif data['m3'][res]['chi2_test'] and data['f-m1m3'][res]['ftest']: 
123                                  self.mf.data.results[res] = self.fill_results(data['m3'][res], model='3') 
124                                  if self.mf.debug: 
125                                          self.mf.log.write('%-12s' % '[Model 3]') 
126   
127                          # Large chi squared test for model 1. 
128                          elif data['m1'][res]['large_chi2'] == 0: 
129                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='1') 
130                                  if self.mf.debug: 
131                                          self.mf.log.write('%-12s' % '[Model 1*]') 
132   
133                          # Test if both model 4 and 5 fit!!! (Should not occur) 
134                          elif data['m4'][res]['zero_chi2'] and data['m5'][res]['zero_chi2']: 
135                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='4+5') 
136                                  if self.mf.debug: 
137                                          self.mf.log.write('%-12s' % '[Model 4 and 5]') 
138   
139                          # Model 4 test. 
140                          elif data['m4'][res]['zero_chi2']: 
141                                  self.mf.data.results[res] = self.fill_results(data['m4'][res], model='4') 
142                                  if self.mf.debug: 
143                                          self.mf.log.write('%-12s' % '[Model 4]') 
144   
145                          # Model 5 test. 
146                          elif data['m5'][res]['zero_chi2']: 
147                                  self.mf.data.results[res] = self.fill_results(data['m5'][res], model='5') 
148                                  if self.mf.debug: 
149                                          self.mf.log.write('%-12s' % '[Model 5]') 
150   
151                          # No model fits! 
152                          else: 
153                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='0') 
154                                  if self.mf.debug: 
155                                          self.mf.log.write('%-12s' % '[Model 0]') 
156   
157   
158 -        def model_selection_extended(self): 
159                  """Palmer's model selection (extended). 
160   
161                  Palmer's model selection, but with additional chi-squared and F-tests for models 4 and 5 
162                  because the number of degrees of freedom for these models are greater than 0.  See the code 
163                  below for details of these changes. 
164                  """ 
165   
166                  data = self.mf.data.data 
167   
168                  if self.mf.debug: 
169                          self.mf.log.write("\n\n<<< Palmer's model selection (extended) >>>\n\n") 
170   
171                  for res in range(len(self.mf.data.relax_data[0])): 
172                          self.mf.data.results.append({}) 
173   
174                          if self.mf.debug: 
175                                  self.mf.log.write('%-22s\n' % ( "Checking res " + data['m1'][res]['res_num'] )) 
176   
177                          # Model 1 test. 
178                          if data['m1'][res]['chi2_test']: 
179                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='1') 
180                                  if self.mf.debug: 
181                                          self.mf.log.write('%-12s' % '[Model 1]') 
182   
183                          # Test if both model 2 and 3 fit!!! (Should not occur) 
184                          elif data['m2'][res]['chi2_test'] and data['f-m1m2'][res]['ftest'] \ 
185                                  and data['m3'][res]['chi2_test'] and data['f-m1m3'][res]['ftest']: 
186                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='2+3') 
187                                  if self.mf.debug: 
188                                          self.mf.log.write('%-12s' % '[Model 2 and 3]') 
189   
190                          # Model 2 test. 
191                          elif data['m2'][res]['chi2_test'] and data['f-m1m2'][res]['ftest']: 
192                                  self.mf.data.results[res] = self.fill_results(data['m2'][res], model='2') 
193                                  if self.mf.debug: 
194                                          self.mf.log.write('%-12s' % '[Model 2]') 
195   
196                          # Model 3 test. 
197                          elif data['m3'][res]['chi2_test'] and data['f-m1m3'][res]['ftest']: 
198                                  self.mf.data.results[res] = self.fill_results(data['m3'][res], model='3') 
199                                  if self.mf.debug: 
200                                          self.mf.log.write('%-12s' % '[Model 3]') 
201   
202                          # Large chi squared test for model 1. 
203                          elif data['m1'][res]['large_chi2'] == 0: 
204                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='1') 
205                                  if self.mf.debug: 
206                                          self.mf.log.write('%-12s' % '[Model 1*]') 
207   
208                          # Test if both model 4 and 5 fit!!! (Should not occur) 
209                          elif data['m4'][res]['chi2_test'] and ( data['f-m2m4'][res]['ftest'] or data['f-m3m4'][res]['ftest'] ) \ 
210                                  and data['m5'][res]['chi2_test'] and data['f-m2m5'][res]['ftest']: 
211                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='4+5') 
212                                  if self.mf.debug: 
213                                          self.mf.log.write('%-12s' % '[Model 4 and 5]') 
214   
215                          # Model 4 test. 
216                          elif data['m4'][res]['chi2_test'] and ( data['f-m2m4'][res]['ftest'] or data['f-m3m4'][res]['ftest'] ): 
217                                  self.mf.data.results[res] = self.fill_results(data['m4'][res], model='4') 
218                                  if self.mf.debug: 
219                                          self.mf.log.write('%-12s' % '[Model 4]') 
220   
221                          # Model 5 test. 
222                          elif data['m5'][res]['chi2_test'] and data['f-m2m5'][res]['ftest']: 
223                                  self.mf.data.results[res] = self.fill_results(data['m5'][res], model='5') 
224                                  if self.mf.debug: 
225                                          self.mf.log.write('%-12s' % '[Model 5]') 
226   
227                          # No model fits! 
228                          else: 
229                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='0') 
230                                  if self.mf.debug: 
231                                          self.mf.log.write('%-12s' % '[Model 0]') 
232   
233   
234 -        def print_data(self): 
235                  "Print all the data into the 'data_all' file." 
236   
237                  file = open('data_all', 'w') 
238   
239                  sys.stdout.write("[") 
240                  for res in range(len(self.mf.data.results)): 
241                          sys.stdout.write("-") 
242                          file.write("\n\n<<< Residue " + self.mf.data.results[res]['res_num']) 
243                          file.write(", Model " + self.mf.data.results[res]['model'] + " >>>\n") 
244                          file.write('%-20s' % '') 
245                          file.write('%-19s' % 'Model 1') 
246                          file.write('%-19s' % 'Model 2') 
247                          file.write('%-19s' % 'Model 3') 
248                          file.write('%-19s' % 'Model 4') 
249                          file.write('%-19s' % 'Model 5') 
250   
251                          # S2. 
252                          file.write('\n%-20s' % 'S2') 
253                          for run in self.mf.data.runs: 
254                                  if match('^m', run): 
255                                          file.write('%9.3f' % self.mf.data.data[run][res]['s2']) 
256                                          file.write('%1s' % '±') 
257                                          file.write('%-9.3f' % self.mf.data.data[run][res]['s2_err']) 
258   
259                          # S2f. 
260                          file.write('\n%-20s' % 'S2f') 
261                          for run in self.mf.data.runs: 
262                                  if match('^m', run): 
263                                          file.write('%9.3f' % self.mf.data.data[run][res]['s2f']) 
264                                          file.write('%1s' % '±') 
265                                          file.write('%-9.3f' % self.mf.data.data[run][res]['s2f_err']) 
266   
267                          # S2s. 
268                          file.write('\n%-20s' % 'S2s') 
269                          for run in self.mf.data.runs: 
270                                  if match('^m', run): 
271                                          file.write('%9.3f' % self.mf.data.data[run][res]['s2s']) 
272                                          file.write('%1s' % '±') 
273                                          file.write('%-9.3f' % self.mf.data.data[run][res]['s2s_err']) 
274   
275                          # te. 
276                          file.write('\n%-20s' % 'te') 
277                          for run in self.mf.data.runs: 
278                                  if match('^m', run): 
279                                          file.write('%9.2f' % self.mf.data.data[run][res]['te']) 
280                                          file.write('%1s' % '±') 
281                                          file.write('%-9.2f' % self.mf.data.data[run][res]['te_err']) 
282   
283                          # Rex. 
284                          file.write('\n%-20s' % 'Rex') 
285                          for run in self.mf.data.runs: 
286                                  if match('^m', run): 
287                                          file.write('%9.3f' % self.mf.data.data[run][res]['rex']) 
288                                          file.write('%1s' % '±') 
289                                          file.write('%-9.3f' % self.mf.data.data[run][res]['rex_err']) 
290   
291                          # Chi2. 
292                          file.write('\n%-20s' % 'Chi2') 
293                          for run in self.mf.data.runs: 
294                                  if match('^m', run): 
295                                          file.write('%-19.3f' % self.mf.data.data[run][res]['chi2']) 
296   
297                          # Chi2 limit. 
298                          file.write('\n%-20s' % 'Chi2 limit') 
299                          for run in self.mf.data.runs: 
300                                  if match('^m', run): 
301                                          file.write('%-19.3f' % self.mf.data.data[run][res]['chi2_lim']) 
302   
303                          # Chi2 test. 
304                          file.write('\n%-20s' % 'Chi2 test') 
305                          for run in self.mf.data.runs: 
306                                  if match('^m', run): 
307                                          file.write('%-19s' % self.mf.data.data[run][res]['chi2_test']) 
308   
309                          # large Chi2. 
310                          file.write('\n%-20s' % 'large Chi2') 
311                          for run in self.mf.data.runs: 
312                                  if match('^m', run): 
313                                          file.write('%-19s' % self.mf.data.data[run][res]['large_chi2']) 
314   
315                          # zero Chi2. 
316                          file.write('\n%-20s' % 'zero Chi2') 
317                          for run in self.mf.data.runs: 
318                                  if match('^m', run): 
319                                          file.write('%-19s' % self.mf.data.data[run][res]['zero_chi2']) 
320   
321                  file.write('\n') 
322                  sys.stdout.write("]\n") 
323                  file.close() 
324