Source Code for Module asymptotic

  1  # A method based on asymptotic model selection criteria. 
  2  # 
  3  # The following asymptotic methods are supported: 
  4  #       AIC - Akaike Information Criteria 
  5  #       BIC - Schwartz Criteria 
  6  # 
  7  # The program is divided into the following stages: 
  8  #       Stage 1:  Creation of the files for the model-free calculations for models 1 to 5.  Monte Carlo 
  9  #               simulations are not used on these initial runs, because the errors are not needed (should 
 10  #               speed up analysis considerably). 
 11  #       Stage 2:  Model selection and the creation of the final run.  Monte Carlo simulations are used to 
 12  #               find errors.  This stage has the option of optimizing the diffusion tensor along with the 
 13  #               model-free parameters. 
 14  #       Stage 3:  Extraction of the data. 
 15   
 16  import sys 
 17  from math import log, pi 
 18  from re import match 
 19   
 20  from common_ops import common_operations 
 21   
 22   
 23 -class asymptotic(common_operations): 
 24 -        def __init__(self, mf): 
 25                  "Model-free analysis based on asymptotic model selection methods." 
 26   
 27                  self.mf = mf 
 28   
 29                  print "Model-free analysis based on " + self.mf.data.usr_param.method + " model selection." 
 30                  self.initialize() 
 31                  self.mf.data.runs = ['m1', 'm2', 'm3', 'm4', 'm5'] 
 32                  self.goto_stage() 
 33   
 34   
 35 -        def calc_crit(self, res, n, k, chisq): 
 36                  "Calculate the criteria" 
 37   
 38                  sum_ln_err = 0 
 39                  for i in range(len(self.mf.data.relax_data)): 
 40                          if self.mf.data.relax_data[i][res][3] == 0: 
 41                                  ln_err = -99 
 42                          else: 
 43                                  ln_err = log(float(self.mf.data.relax_data[i][res][3])) 
 44                          sum_ln_err = sum_ln_err + ln_err 
 45   
 46                  if match('^AIC$', self.mf.data.usr_param.method): 
 47                          aic = n*log(2*pi) + sum_ln_err + chisq + 2*k 
 48                          aic = aic / (2*n) 
 49                          return aic 
 50   
 51                  elif match('^AICc$', self.mf.data.usr_param.method): 
 52                          aicc = n*log(2*pi) + sum_ln_err + chisq + 2*k + 2*k*(k+1)/(n-k-1) 
 53                          aicc = aicc / (2*n) 
 54                          return aicc 
 55   
 56                  elif match('^BIC$', self.mf.data.usr_param.method): 
 57                          bic = n*log(2*pi) + sum_ln_err + chisq + k*log(n) 
 58                          bic = bic / ( 2 * n ) 
 59                          return bic 
 60   
 61   
 62 -        def initial_runs(self): 
 63                  "Creation of the files for the Modelfree calculations for models 1 to 5." 
 64                   
 65                  for run in self.mf.data.runs: 
 66                          print "Creating input files for model " + run 
 67                          self.mf.log.write("\n\n<<< Model " + run + " >>>\n\n") 
 68                          self.mf.file_ops.mkdir(dir=run) 
 69                          self.mf.file_ops.open_mf_files(dir=run) 
 70                          self.set_run_flags(run) 
 71                          self.log_params('M1', self.mf.data.usr_param.md1) 
 72                          self.log_params('M2', self.mf.data.usr_param.md2) 
 73                          self.create_mfin(sims='n', sim_type='pred') 
 74                          self.create_run(dir=run) 
 75                          for res in range(len(self.mf.data.relax_data[0])): 
 76                                  # Mfdata. 
 77                                  self.create_mfdata(res) 
 78                                  # Mfmodel. 
 79                                  self.create_mfmodel(res, self.mf.data.usr_param.md1, type='M1') 
 80                                  # Mfpar. 
 81                                  self.create_mfpar(res) 
 82                          self.mf.file_ops.close_mf_files(dir=run) 
 83   
 84   
 85 -        def model_selection(self): 
 86                  "Model selection." 
 87                   
 88                  data = self.mf.data.data 
 89   
 90                  self.mf.log.write("\n\n<<< " + self.mf.data.usr_param.method + " model selection >>>") 
 91                  for res in range(len(self.mf.data.relax_data[0])): 
 92                          self.mf.data.results.append({}) 
 93                          self.mf.log.write('\n%-22s' % ( "   Checking res " + data['m1'][res]['res_num'] )) 
 94   
 95                          n = self.mf.data.num_data_sets 
 96   
 97                          data['m1'][res]['crit'] = self.calc_crit(res, n, k=1, chisq=data['m1'][res]['sse']) 
 98                          data['m2'][res]['crit'] = self.calc_crit(res, n, k=2, chisq=data['m2'][res]['sse']) 
 99                          data['m3'][res]['crit'] = self.calc_crit(res, n, k=2, chisq=data['m3'][res]['sse']) 
100                          data['m4'][res]['crit'] = self.calc_crit(res, n, k=3, chisq=data['m4'][res]['sse']) 
101                          data['m5'][res]['crit'] = self.calc_crit(res, n, k=3, chisq=data['m5'][res]['sse']) 
102   
103                          # Select model. 
104                          min = 'm1' 
105                          for run in self.mf.data.runs: 
106                                  if data[run][res]['crit'] < data[min][res]['crit']: 
107                                          min = run 
108                          self.mf.data.results[res] = self.fill_results(data[min][res], model=min[1]) 
109   
110                          self.mf.log.write("\n\t" + self.mf.data.usr_param.method + " (m1): " + `data['m1'][res]['crit']` + "\n") 
111                          self.mf.log.write("\n\t" + self.mf.data.usr_param.method + " (m2): " + `data['m2'][res]['crit']` + "\n") 
112                          self.mf.log.write("\n\t" + self.mf.data.usr_param.method + " (m3): " + `data['m3'][res]['crit']` + "\n") 
113                          self.mf.log.write("\n\t" + self.mf.data.usr_param.method + " (m4): " + `data['m4'][res]['crit']` + "\n") 
114                          self.mf.log.write("\n\t" + self.mf.data.usr_param.method + " (m5): " + `data['m5'][res]['crit']` + "\n") 
115                          self.mf.log.write("\tThe selected model is: " + min + "\n\n") 
116   
117   
118 -        def stage2(self): 
119                  self.mf.file_ops.mkdir('grace') 
120                   
121                  print "\n[ Model-free data extraction ]\n" 
122                  for run in self.mf.data.runs: 
123                          mfout = self.mf.file_ops.read_file(run + '/mfout') 
124                          mfout_lines = mfout.readlines() 
125                          mfout.close() 
126                          print "Extracting model-free data from " + run + "/mfout." 
127                          num_res = len(self.mf.data.relax_data[0]) 
128                          self.mf.data.data[run] = self.mf.star.extract(mfout_lines, num_res) 
129   
130                  print "\n[ " + self.mf.data.usr_param.method + " model selection ]\n" 
131                  self.model_selection() 
132   
133                  print "\n[ Printing results ]\n" 
134                  self.print_results() 
135   
136                  print "\n[ Placing data structures into \"data_all\" ]\n" 
137                  self.print_data() 
138   
139                  print "\n[ Grace file creation ]\n" 
140                  self.grace('grace/S2.agr', 'S2', subtitle="After model selection, unoptimized") 
141                  self.grace('grace/S2f.agr', 'S2f', subtitle="After model selection, unoptimized") 
142                  self.grace('grace/S2s.agr', 'S2s', subtitle="After model selection, unoptimized") 
143                  self.grace('grace/te.agr', 'te', subtitle="After model selection, unoptimized") 
144                  self.grace('grace/Rex.agr', 'Rex', subtitle="After model selection, unoptimized") 
145                  self.grace('grace/SSE.agr', 'SSE', subtitle="After model selection, unoptimized") 
146