Source Code for Module model_selection.asymptotic

  1  # A method based on asymptotic model selection criteria. 
  2  # 
  3  # The following asymptotic methods are supported: 
  4  #       AIC - Akaike Information Criteria 
  5  #       AICc - Akaike Information Criteria corrected for small sample sizes 
  6  #       BIC - Schwartz Criteria 
  7  # 
  8  # The program is divided into the following stages: 
  9  #       Stage 1:  Creation of the files for the model-free calculations for models 1 to 5.  Monte Carlo 
 10  #               simulations are not used on these initial runs, because the errors are not needed (should 
 11  #               speed up analysis considerably). 
 12  #       Stage 2:  Model selection and the creation of the final run.  Monte Carlo simulations are used to 
 13  #               find errors.  This stage has the option of optimizing the diffusion tensor along with the 
 14  #               model-free parameters. 
 15  #       Stage 3:  Extraction of the data. 
 16   
 17  import sys 
 18  from math import log, pi 
 19  from re import match 
 20   
 21  from common_ops import common_operations 
 22   
 23   
 24 -class asymptotic(common_operations): 
 25 -        def __init__(self, mf): 
 26                  "Model-free analysis based on asymptotic model selection methods." 
 27   
 28                  self.mf = mf 
 29   
 30   
 31 -        def run(self): 
 32                  "Model selection." 
 33   
 34                  data = self.mf.data.data 
 35                  n = float(self.mf.data.num_data_sets) 
 36   
 37                  if self.mf.debug: 
 38                          self.mf.log.write("\n\n<<< " + self.mf.usr_param.method + " model selection >>>\n\n") 
 39   
 40                  for res in range(len(self.mf.data.relax_data[0])): 
 41                          self.mf.data.results.append({}) 
 42   
 43                          if self.mf.debug: 
 44                                  self.mf.log.write('%-22s\n' % ( "Checking res " + data['m1'][res]['res_num'] )) 
 45   
 46                          err = [] 
 47                          for set in range(len(self.mf.data.relax_data)): 
 48                                  err.append(float(self.mf.data.relax_data[set][res][3])) 
 49   
 50                          for model in self.mf.data.runs: 
 51                                  chi2 = data[model][res]['chi2'] 
 52                                  crit = chi2 / (2.0 * n) 
 53   
 54                                  if match('m1', model): 
 55                                          k = 1.0 
 56                                  elif match('m2', model) or match('m3', model): 
 57                                          k = 2.0 
 58                                  elif match('m4', model) or match('m5', model): 
 59                                          k = 3.0 
 60   
 61                                  if match('^AIC$', self.mf.usr_param.method): 
 62                                          crit = crit + k / n 
 63   
 64                                  elif match('^AICc$', self.mf.usr_param.method): 
 65                                          if n - k == 1: 
 66                                                  crit = 1e99 
 67                                          else: 
 68                                                  crit = crit + k/n + k*(k + 1.0)/((n - k - 1.0) * n) 
 69   
 70                                  elif match('^BIC$', self.mf.usr_param.method): 
 71                                          crit = crit + k*log(n) / (2.0 * n) 
 72   
 73                                  data[model][res]['crit'] = crit 
 74   
 75                          # Select model. 
 76                          min = 'm1' 
 77                          for run in self.mf.data.runs: 
 78                                  if data[run][res]['crit'] < data[min][res]['crit']: 
 79                                          min = run 
 80                          if data[min][res]['crit'] == float('inf'): 
 81                                  self.mf.data.results[res] = self.fill_results(data[min][res], model='0') 
 82                          else: 
 83                                  self.mf.data.results[res] = self.fill_results(data[min][res], model=min[1]) 
 84   
 85                          if self.mf.debug: 
 86                                  self.mf.log.write(self.mf.usr_param.method + " (m1): " + `data['m1'][res]['crit']` + "\n") 
 87                                  self.mf.log.write(self.mf.usr_param.method + " (m2): " + `data['m2'][res]['crit']` + "\n") 
 88                                  self.mf.log.write(self.mf.usr_param.method + " (m3): " + `data['m3'][res]['crit']` + "\n") 
 89                                  self.mf.log.write(self.mf.usr_param.method + " (m4): " + `data['m4'][res]['crit']` + "\n") 
 90                                  self.mf.log.write(self.mf.usr_param.method + " (m5): " + `data['m5'][res]['crit']` + "\n") 
 91                                  self.mf.log.write("The selected model is: " + min + "\n\n") 
 92   
 93   
 94 -        def print_data(self): 
 95                  "Print all the data into the 'data_all' file." 
 96   
 97                  file = open('data_all', 'w') 
 98                  file_crit = open('crit', 'w') 
 99   
100                  sys.stdout.write("[") 
101                  for res in range(len(self.mf.data.results)): 
102                          sys.stdout.write("-") 
103                          file.write("\n\n<<< Residue " + self.mf.data.results[res]['res_num']) 
104                          file.write(", Model " + self.mf.data.results[res]['model'] + " >>>\n") 
105                          file.write('%-20s' % '') 
106                          file.write('%-19s' % 'Model 1') 
107                          file.write('%-19s' % 'Model 2') 
108                          file.write('%-19s' % 'Model 3') 
109                          file.write('%-19s' % 'Model 4') 
110                          file.write('%-19s' % 'Model 5') 
111   
112                          file_crit.write('%-6s' % self.mf.data.results[res]['res_num']) 
113                          file_crit.write('%-6s' % self.mf.data.results[res]['model']) 
114   
115                          # S2. 
116                          file.write('\n%-20s' % 'S2') 
117                          for run in self.mf.data.runs: 
118                                  file.write('%9.3f' % self.mf.data.data[run][res]['s2']) 
119                                  file.write('%1s' % '±') 
120                                  file.write('%-9.3f' % self.mf.data.data[run][res]['s2_err']) 
121   
122                          # S2f. 
123                          file.write('\n%-20s' % 'S2f') 
124                          for run in self.mf.data.runs: 
125                                  file.write('%9.3f' % self.mf.data.data[run][res]['s2f']) 
126                                  file.write('%1s' % '±') 
127                                  file.write('%-9.3f' % self.mf.data.data[run][res]['s2f_err']) 
128   
129                          # S2s. 
130                          file.write('\n%-20s' % 'S2s') 
131                          for run in self.mf.data.runs: 
132                                  file.write('%9.3f' % self.mf.data.data[run][res]['s2s']) 
133                                  file.write('%1s' % '±') 
134                                  file.write('%-9.3f' % self.mf.data.data[run][res]['s2s_err']) 
135   
136                          # te. 
137                          file.write('\n%-20s' % 'te') 
138                          for run in self.mf.data.runs: 
139                                  file.write('%9.2f' % self.mf.data.data[run][res]['te']) 
140                                  file.write('%1s' % '±') 
141                                  file.write('%-9.2f' % self.mf.data.data[run][res]['te_err']) 
142   
143                          # Rex. 
144                          file.write('\n%-20s' % 'Rex') 
145                          for run in self.mf.data.runs: 
146                                  file.write('%9.3f' % self.mf.data.data[run][res]['rex']) 
147                                  file.write('%1s' % '±') 
148                                  file.write('%-9.3f' % self.mf.data.data[run][res]['rex_err']) 
149   
150                          # Chi2. 
151                          file.write('\n%-20s' % 'Chi2') 
152                          for run in self.mf.data.runs: 
153                                  file.write('%-19.3f' % self.mf.data.data[run][res]['chi2']) 
154   
155                          # Model selection criteria. 
156                          file.write('\n%-20s' % self.mf.usr_param.method) 
157                          for run in self.mf.data.runs: 
158                                  file.write('%-19.6f' % self.mf.data.data[run][res]['crit']) 
159   
160                                  file_crit.write('%-25s' % `self.mf.data.data[run][res]['crit']`) 
161                          file_crit.write('\n') 
162   
163                  file.write('\n') 
164                  sys.stdout.write("]\n") 
165                  file.close() 
166