Source Code for Module 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                  self.mf = mf 
 23   
 24                  print "Palmer's method for model-free analysis. (Mandel et al., 1995)" 
 25                  self.initialize() 
 26                  self.mf.data.runs = ['m1', 'm2', 'm3', 'm4', 'm5', 'f-m1m2', 'f-m1m3'] 
 27                  if self.mf.data.num_data_sets > 3: 
 28                          print "\nThe number of data sets is greater than 3.  Setting up additional F-tests." 
 29                          self.mf.data.runs.append('f-m2m4') 
 30                          self.mf.data.runs.append('f-m2m5') 
 31                          self.mf.data.runs.append('f-m3m4') 
 32                  self.goto_stage() 
 33   
 34   
 35 -        def initial_runs(self): 
 36                  "Creation of the files for the Modelfree calculations for models 1 to 5 and the F-tests." 
 37                   
 38                  for run in self.mf.data.runs: 
 39                          if match('^m', run): 
 40                                  print "Creating input files for model " + run 
 41                                  self.mf.log.write("\n\n<<< Model " + run + " >>>\n\n") 
 42                          elif match('^f', run): 
 43                                  print "Creating input files for the F-test " + run 
 44                                  self.mf.log.write("\n\n<<< F-test " + run + " >>>\n\n") 
 45                          else: 
 46                                  print "The run '" + run + "'does not start with an m or f, quitting script!\n\n" 
 47                                  sys.exit() 
 48                          self.mf.file_ops.mkdir(dir=run) 
 49                          self.mf.file_ops.open_mf_files(dir=run) 
 50                          self.set_run_flags(run) 
 51                          self.log_params('M1', self.mf.data.usr_param.md1) 
 52                          self.log_params('M2', self.mf.data.usr_param.md2) 
 53                          if match('^m', run): 
 54                                  self.create_mfin(sims='y', sim_type='pred') 
 55                          elif match('^f', run): 
 56                                  self.create_mfin(sel='ftest', sims='y', sim_type='pred') 
 57                          self.create_run(dir=run) 
 58                          for res in range(len(self.mf.data.relax_data[0])): 
 59                                  # Mfdata. 
 60                                  self.create_mfdata(res) 
 61                                  # Mfmodel. 
 62                                  self.create_mfmodel(res, self.mf.data.usr_param.md1, type='M1') 
 63                                  if match('^f', run): 
 64                                          self.create_mfmodel(res, self.mf.data.usr_param.md2, type='M2') 
 65                                  # Mfpar. 
 66                                  self.create_mfpar(res) 
 67                          self.mf.file_ops.close_mf_files(dir=run) 
 68   
 69   
 70 -        def model_selection(self): 
 71                  "Palmer's model selection." 
 72   
 73                  data = self.mf.data.data 
 74   
 75                  self.mf.log.write("\n\n<<< Palmer's model selection >>>") 
 76                  for res in range(len(self.mf.data.relax_data[0])): 
 77                          self.mf.data.results.append({}) 
 78                          self.mf.log.write('\n%-22s' % ( "   Checking res " + data['m1'][res]['res_num'] )) 
 79   
 80                          # Model 1 test. 
 81                          if data['m1'][res]['sse_test'] == 1: 
 82                                  self.mf.log.write('%-12s' % '[Model 1]') 
 83                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='1') 
 84   
 85                          # Test if both model 2 and 3 fit!!! (Should not occur) 
 86                          elif data['m2'][res]['sse_test'] == 1 and data['f-m1m2'][res]['ftest'] == 1 \ 
 87                                  and data['m3'][res]['sse_test'] == 1 and data['f-m1m3'][res]['ftest'] == 1: 
 88                                  self.mf.log.write('%-12s' % '[Model 2 and 3]') 
 89                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='2+3') 
 90   
 91                          # Model 2 test. 
 92                          elif data['m2'][res]['sse_test'] == 1 and data['f-m1m2'][res]['ftest'] == 1: 
 93                                  self.mf.log.write('%-12s' % '[Model 2]') 
 94                                  self.mf.data.results[res] = self.fill_results(data['m2'][res], model='2') 
 95   
 96                          # Model 3 test. 
 97                          elif data['m3'][res]['sse_test'] == 1 and data['f-m1m3'][res]['ftest'] == 1: 
 98                                  self.mf.log.write('%-12s' % '[Model 3]') 
 99                                  self.mf.data.results[res] = self.fill_results(data['m3'][res], model='3') 
100   
101                          # Large SSE test for model 1. 
102                          elif data['m1'][res]['large_sse'] == 0: 
103                                  self.mf.log.write('%-12s' % '[Model 1*]') 
104                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='1') 
105   
106                          # Test if both model 4 and 5 fit!!! (Should not occur) 
107                          elif data['m4'][res]['zero_sse'] == 1 and data['m5'][res]['zero_sse'] == 1: 
108                                  self.mf.log.write('%-12s' % '[Model 4 and 5]') 
109                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='4+5') 
110                                   
111                          # Model 4 test. 
112                          elif data['m4'][res]['zero_sse'] == 1: 
113                                  self.mf.log.write('%-12s' % '[Model 4]') 
114                                  self.mf.data.results[res] = self.fill_results(data['m4'][res], model='4') 
115   
116                          # Model 5 test. 
117                          elif data['m5'][res]['zero_sse'] == 1: 
118                                  self.mf.log.write('%-12s' % '[Model 5]') 
119                                  self.mf.data.results[res] = self.fill_results(data['m5'][res], model='5') 
120   
121                          # No model fits! 
122                          else: 
123                                  self.mf.log.write('%-12s' % '[Model 0]') 
124                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='0') 
125   
126   
127 -        def model_selection_extended(self): 
128                  """Palmer's model selection (extended). 
129                   
130                  Palmer's model selection, but with additional chi-squared and F-tests for models 4 and 5 
131                  because the number of degrees of freedom for these models are greater than 0.  See the code 
132                  below for details of these changes. 
133                  """ 
134   
135                  data = self.mf.data.data 
136   
137                  self.mf.log.write("\n\n<<< Palmer's model selection (extended) >>>") 
138                  for res in range(len(self.mf.data.relax_data[0])): 
139                          self.mf.data.results.append({}) 
140                          self.mf.log.write('\n%-22s' % ( "   Checking res " + data['m1'][res]['res_num'] )) 
141   
142                          # Model 1 test. 
143                          if data['m1'][res]['sse_test'] == 1: 
144                                  self.mf.log.write('%-12s' % '[Model 1]') 
145                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='1') 
146   
147                          # Test if both model 2 and 3 fit!!! (Should not occur) 
148                          elif data['m2'][res]['sse_test'] == 1 and data['f-m1m2'][res]['ftest'] == 1 \ 
149                                  and data['m3'][res]['sse_test'] == 1 and data['f-m1m3'][res]['ftest'] == 1: 
150                                  self.mf.log.write('%-12s' % '[Model 2 and 3]') 
151                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='2+3') 
152   
153                          # Model 2 test. 
154                          elif data['m2'][res]['sse_test'] == 1 and data['f-m1m2'][res]['ftest'] == 1: 
155                                  self.mf.log.write('%-12s' % '[Model 2]') 
156                                  self.mf.data.results[res] = self.fill_results(data['m2'][res], model='2') 
157   
158                          # Model 3 test. 
159                          elif data['m3'][res]['sse_test'] == 1 and data['f-m1m3'][res]['ftest'] == 1: 
160                                  self.mf.log.write('%-12s' % '[Model 3]') 
161                                  self.mf.data.results[res] = self.fill_results(data['m3'][res], model='3') 
162   
163                          # Large SSE test for model 1. 
164                          elif data['m1'][res]['large_sse'] == 0: 
165                                  self.mf.log.write('%-12s' % '[Model 1*]') 
166                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='1') 
167   
168                          # Test if both model 4 and 5 fit!!! (Should not occur) 
169                          elif data['m4'][res]['sse_test'] == 1 and ( data['f-m2m4'][res]['ftest'] == 1 or data['f-m3m4'][res]['ftest'] == 1 ) \ 
170                                  and data['m5'][res]['sse_test'] == 1 and data['f-m2m5'][res]['ftest'] == 1: 
171                                  self.mf.log.write('%-12s' % '[Model 4 and 5]') 
172                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='4+5') 
173                                   
174                          # Model 4 test. 
175                          elif data['m4'][res]['sse_test'] == 1 and ( data['f-m2m4'][res]['ftest'] == 1 or data['f-m3m4'][res]['ftest'] == 1 ): 
176                                  self.mf.log.write('%-12s' % '[Model 4]') 
177                                  self.mf.data.results[res] = self.fill_results(data['m4'][res], model='4') 
178                                   
179                          # Model 5 test. 
180                          elif data['m5'][res]['sse_test'] == 1 and data['f-m2m5'][res]['ftest'] == 1: 
181                                  self.mf.log.write('%-12s' % '[Model 5]') 
182                                  self.mf.data.results[res] = self.fill_results(data['m5'][res], model='5') 
183   
184                          # No model fits! 
185                          else: 
186                                  self.mf.log.write('%-12s' % '[Model 0]') 
187                                  self.mf.data.results[res] = self.fill_results(data['m1'][res], model='0') 
188   
189   
190 -        def stage2(self): 
191                  self.mf.file_ops.mkdir('grace') 
192   
193                  print "\n[ Model-free data extraction ]\n" 
194                  for run in self.mf.data.runs: 
195                          mfout = self.mf.file_ops.read_file(run + '/mfout') 
196                          mfout_lines = mfout.readlines() 
197                          mfout.close() 
198                          print "Extracting model-free data from " + run + "/mfout." 
199                          num_res = len(self.mf.data.relax_data[0]) 
200                          if match('^m', run): 
201                                  self.mf.data.data[run] = self.mf.star.extract(mfout_lines, num_res, self.mf.data.usr_param.sse_lim, self.mf.data.usr_param.ftest_lim, float(self.mf.data.usr_param.large_sse), ftest='n') 
202                          if match('^f', run): 
203                                  self.mf.data.data[run] = self.mf.star.extract(mfout_lines, num_res, self.mf.data.usr_param.sse_lim, self.mf.data.usr_param.ftest_lim, float(self.mf.data.usr_param.large_sse), ftest='y') 
204   
205                  print "\n[ Palmer's model selection ]\n" 
206                  if self.mf.data.num_data_sets > 3: 
207                          # ie degrees of freedom > 0 in all models. 
208                          self.mf.log.write("Extended model selection.\n\n") 
209                          print "The number of data sets is greater than 3." 
210                          print "\tRunning Palmer's model selection, with additional chi-squared and F-tests" 
211                          print "\tfor models 4 and 5 (the degrees of freedom for these models are greater than 0).\n" 
212                          self.model_selection_extended() 
213                  else: 
214                          self.mf.log.write("Normal model selection.\n\n") 
215                          print "The number of data sets is equal to 3." 
216                          print "\tRunning Palmer's model selection.\n" 
217                          self.model_selection() 
218   
219                  print "\n[ Printing results ]\n" 
220                  self.print_results() 
221   
222                  print "\n[ Placing data structures into \"data_all\" ]\n" 
223                  self.print_data(ftests='y') 
224   
225                  print "\n[ Grace file creation ]\n" 
226                  self.grace('grace/S2.agr', 'S2', subtitle="After model selection, unoptimized") 
227                  self.grace('grace/S2f.agr', 'S2f', subtitle="After model selection, unoptimized") 
228                  self.grace('grace/S2s.agr', 'S2s', subtitle="After model selection, unoptimized") 
229                  self.grace('grace/te.agr', 'te', subtitle="After model selection, unoptimized") 
230                  self.grace('grace/Rex.agr', 'Rex', subtitle="After model selection, unoptimized") 
231                  self.grace('grace/SSE.agr', 'SSE', subtitle="After model selection, unoptimized") 
232