Source Code for Module auto_analyses.dauvergne_protocol

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2004-2012 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the program relax.                                     # 
  6  #                                                                             # 
  7  # relax is free software; you can redistribute it and/or modify               # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation; either version 2 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # relax is distributed in the hope that it will be useful,                    # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with relax; if not, write to the Free Software                        # 
 19  # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA   # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23  # Python module imports. 
 24  from os import F_OK, access, getcwd, listdir, sep 
 25  from re import search 
 26  from string import lower 
 27  from time import sleep 
 28   
 29  # relax module imports. 
 30  from doc_builder import LIST, PARAGRAPH, SECTION, SUBSECTION, TITLE, to_docstring 
 31  from float import floatAsByteArray 
 32  from info import Info_box; info = Info_box() 
 33  from generic_fns.mol_res_spin import exists_mol_res_spin_data, generate_spin_id, spin_index_loop, spin_loop 
 34  from generic_fns.pipes import cdp_name, get_pipe, has_pipe, pipe_names, switch 
 35  from generic_fns import selection 
 36  from prompt.interpreter import Interpreter 
 37  from relax_errors import RelaxError, RelaxNoSequenceError, RelaxNoValueError 
 38  from relax_io import DummyFileObject 
 39  from status import Status; status = Status() 
 40   
 41   
 42  doc = [ 
 43          [TITLE, "Automatic analysis for black-box model-free results."], 
 44          [PARAGRAPH, "The dauvergne_protocol auto-analysis is designed for those who appreciate black-boxes or those who appreciate complex code.  Importantly, data at multiple magnetic field strengths is essential for this analysis.  If you would like to change how model-free analysis is performed, the code in the file auto_analyses/dauvergne_protocol.py in the base relax directory can be copied and modified as needed and used with the relax script interface.  This file is simply a complex relax script.  For a description of object-oriented coding in python using classes, functions/methods, self, etc., please see the python tutorial."], 
 45   
 46          [SECTION, "References"], 
 47   
 48          [SUBSECTION, "Auto-analysis primary reference"], 
 49          [PARAGRAPH, "The model-free optimisation methodology herein is that of:"], 
 50          [LIST, info.bib['dAuvergneGooley08b'].cite_short()], 
 51   
 52          [SUBSECTION, "Techniques used in the auto-analysis"], 
 53          [PARAGRAPH, "Other references for features of this dauvergne_protocol auto-analysis include model-free model selection using Akaike's Information Criterion:"], 
 54          [LIST, info.bib['dAuvergneGooley03'].cite_short()], 
 55          [PARAGRAPH, "The elimination of failed model-free models and Monte Carlo simulations:"], 
 56          [LIST, info.bib['dAuvergneGooley06'].cite_short()], 
 57          [PARAGRAPH, "Significant model-free optimisation improvements:"], 
 58          [LIST, info.bib['dAuvergneGooley08a'].cite_short()], 
 59          [PARAGRAPH, "Rather than searching for the lowest chi-squared value, this auto-analysis searches for the model with the lowest AIC criterion.  This complex multi-universe, multi-dimensional problem is formulated, using set theory, as the universal solution:"], 
 60          [LIST, info.bib['dAuvergneGooley07'].cite_short()], 
 61          [PARAGRAPH, "The basic three references for the original and extended model-free theories are:"], 
 62          [LIST, info.bib['LipariSzabo82a'].cite_short()], 
 63          [LIST, info.bib['LipariSzabo82b'].cite_short()], 
 64          [LIST, info.bib['Clore90'].cite_short()], 
 65   
 66          [SECTION, "How to use this auto-analysis"], 
 67          [PARAGRAPH, "The five diffusion models used in this auto-analysis are:"], 
 68          [LIST, "Model I   (MI)   - Local tm."], 
 69          [LIST, "Model II  (MII)  - Sphere."], 
 70          [LIST, "Model III (MIII) - Prolate spheroid."], 
 71          [LIST, "Model IV  (MIV)  - Oblate spheroid."], 
 72          [LIST, "Model V   (MV)   - Ellipsoid."], 
 73          [PARAGRAPH, "If using the script-based user interface (UI), changing the value of the variable diff_model will determine the behaviour of this auto-analysis.  Model I must be optimised prior to any of the other diffusion models, while the Models II to V can be optimised in any order.  To select the various models, set the variable diff_model to the following strings:"], 
 74          [LIST, "MI   - 'local_tm'"], 
 75          [LIST, "MII  - 'sphere'"], 
 76          [LIST, "MIII - 'prolate'"], 
 77          [LIST, "MIV  - 'oblate'"], 
 78          [LIST, "MV   - 'ellipsoid'"], 
 79          [PARAGRAPH, "This approach has the advantage of eliminating the need for an initial estimate of a global diffusion tensor and removing all the problems associated with the initial estimate."], 
 80          [PARAGRAPH, "It is important that the number of parameters in a model does not exceed the number of relaxation data sets for that spin.  If this is the case, the list of models in the mf_models and local_tm_models variables will need to be trimmed."], 
 81   
 82          [SUBSECTION, "Model I - Local tm"], 
 83          [PARAGRAPH, "This will optimise the diffusion model whereby all spin of the molecule have a local tm value, i.e. there is no global diffusion tensor.  This model needs to be optimised prior to optimising any of the other diffusion models.  Each spin is fitted to the multiple model-free models separately, where the parameter tm is included in each model."], 
 84          [PARAGRAPH, "AIC model selection is used to select the models for each spin."], 
 85   
 86          [SUBSECTION, "Model II - Sphere"], 
 87          [PARAGRAPH, "This will optimise the isotropic diffusion model.  Multiple steps are required, an initial optimisation of the diffusion tensor, followed by a repetitive optimisation until convergence of the diffusion tensor.  In the relax script UI each of these steps requires this script to be rerun, unless the conv_loop flag is True.  In the GUI (graphical user interface), the procedure is repeated automatically until convergence.  For the initial optimisation, which will be placed in the directory './sphere/init/', the following steps are used:"], 
 88          [PARAGRAPH, "The model-free models and parameter values for each spin are set to those of diffusion model MI."], 
 89          [PARAGRAPH, "The local tm parameter is removed from the models."], 
 90          [PARAGRAPH, "The model-free parameters are fixed and a global spherical diffusion tensor is minimised."], 
 91          [PARAGRAPH, "For the repetitive optimisation, each minimisation is named from 'round_1' onwards.  The initial 'round_1' optimisation will extract the diffusion tensor from the results file in './sphere/init/', and the results will be placed in the directory './sphere/round_1/'.  Each successive round will take the diffusion tensor from the previous round.  The following steps are used:"], 
 92          [PARAGRAPH, "The global diffusion tensor is fixed and the multiple model-free models are fitted to each spin."], 
 93          [PARAGRAPH, "AIC model selection is used to select the models for each spin."], 
 94          [PARAGRAPH, "All model-free and diffusion parameters are allowed to vary and a global optimisation of all parameters is carried out."], 
 95   
 96          [SUBSECTION, "Model III - Prolate spheroid"], 
 97          [PARAGRAPH, "The methods used are identical to those of diffusion model MII, except that an axially symmetric diffusion tensor with Da >= 0 is used.  The base directory containing all the results is './prolate/'."], 
 98   
 99          [SUBSECTION, "Model IV - Oblate spheroid"], 
100          [PARAGRAPH, "The methods used are identical to those of diffusion model MII, except that an axially symmetric diffusion tensor with Da <= 0 is used.  The base directory containing all the results is './oblate/'."], 
101   
102          [SUBSECTION, "Model V - Ellipsoid"], 
103          [PARAGRAPH, "The methods used are identical to those of diffusion model MII, except that a fully anisotropic diffusion tensor is used (also known as rhombic or asymmetric diffusion).  The base directory is './ellipsoid/'."], 
104   
105          [SUBSECTION, "Final run"], 
106          [PARAGRAPH, "Once all the diffusion models have converged, the final run can be executed.  This is done by setting the variable diff_model to 'final'.  This consists of two steps, diffusion tensor model selection, and Monte Carlo simulations.  Firstly AIC model selection is used to select between the diffusion tensor models.  Monte Carlo simulations are then run solely on this selected diffusion model.  Minimisation of the model is bypassed as it is assumed that the model is already fully optimised (if this is not the case the final run is not yet appropriate)."], 
107          [PARAGRAPH, "The final black-box model-free results will be placed in the file 'final/results'."] 
108  ] 
109   
110   
111  # Build the module docstring. 
112  __doc__ = to_docstring(doc) 
113   
114   
115   
116 -class dAuvergne_protocol: 
117      """The model-free auto-analysis.""" 
118   
119      # Some class variables. 
120      opt_func_tol = 1e-25 
121      opt_max_iterations = int(1e7) 
122   
123 -    def __init__(self, pipe_name=None, results_dir=None, diff_model=None, mf_models=['m0', 'm1', 'm2', 'm3', 'm4', 'm5', 'm6', 'm7', 'm8', 'm9'], local_tm_models=['tm0', 'tm1', 'tm2', 'tm3', 'tm4', 'tm5', 'tm6', 'tm7', 'tm8', 'tm9'], grid_inc=11, diff_tensor_grid_inc={'sphere': 11, 'prolate': 11, 'oblate': 11, 'ellipsoid': 6}, min_algor='newton', mc_sim_num=500, max_iter=None, user_fns=None, conv_loop=True): 
124          """Perform the full model-free analysis protocol of d'Auvergne and Gooley, 2008b. 
125   
126          @keyword pipe_name:             The name of the data pipe containing the sequence info.  This data pipe should have all values set including the CSA value, the bond length, the heteronucleus name and proton name.  It should also have all relaxation data loaded. 
127          @keyword results_dir:           The directory, where files are saved in. 
128          @type results_dir:              str 
129          @keyword diff_model:            The global diffusion model to optimise.  This can be one of 'local_tm', 'sphere', 'oblate', 'prolate', 'ellipsoid', or 'final'.  If all or a subset of these are supplied as a list, then these will be automatically looped over and calculated. 
130          @type diff_model:               str or list of str 
131          @keyword mf_models:             The model-free models. 
132          @type mf_models:                list of str 
133          @keyword local_tm_models:       The model-free models. 
134          @type local_tm_models:          list of str 
135          @keyword grid_inc:              The grid search size (the number of increments per dimension). 
136          @type grid_inc:                 int 
137          @keyword diff_tensor_grid_inc:  A list of grid search sizes for the optimisation of the sphere, prolate spheroid, oblate spheroid, and ellipsoid. 
138          @type diff_tensor_grid_inc:     list of int 
139          @keyword min_algor:             The minimisation algorithm (in most cases this should not be changed). 
140          @type min_algor:                str 
141          @keyword mc_sim_num:            The number of Monte Carlo simulations to be used for error analysis at the end of the analysis. 
142          @type mc_sim_num:               int 
143          @keyword max_iter:              The maximum number of iterations for the global iteration.  Set to None, then the algorithm iterates until convergence. 
144          @type max_iter:                 int or None. 
145          @keyword user_fns:              A dictionary of replacement user functions.  These will overwrite the standard user functions.  The key should be the name of the user function or user function class and the value should be the function or class instance. 
146          @type user_fns:                 dict 
147          @keyword conv_loop:             Automatic looping over all rounds until convergence. 
148          @type conv_loop:                bool 
149          """ 
150   
151          # Execution lock. 
152          status.exec_lock.acquire(pipe_name, mode='auto-analysis') 
153   
154          # Store the args. 
155          self.pipe_name = pipe_name 
156          self.mf_models = mf_models 
157          self.local_tm_models = local_tm_models 
158          self.grid_inc = grid_inc 
159          self.diff_tensor_grid_inc = diff_tensor_grid_inc 
160          self.min_algor = min_algor 
161          self.mc_sim_num = mc_sim_num 
162          self.max_iter = max_iter 
163          self.conv_loop = conv_loop 
164   
165          # The diffusion models. 
166          if isinstance(diff_model, list): 
167              self.diff_model_list = diff_model 
168          else: 
169              self.diff_model_list = [diff_model] 
170   
171          # Project directory (i.e. directory containing the model-free model results and the newly generated files) 
172          if results_dir: 
173              self.results_dir = results_dir + sep 
174          else: 
175              self.results_dir = getcwd() + sep 
176   
177          # Data checks. 
178          self.check_vars() 
179   
180          # Set the data pipe to the current data pipe. 
181          if self.pipe_name != cdp_name(): 
182              switch(self.pipe_name) 
183   
184          # Some info for the status. 
185          self.status_setup() 
186   
187          # Load the interpreter. 
188          self.interpreter = Interpreter(show_script=False, quit=False, raise_relax_error=True) 
189          self.interpreter.populate_self() 
190          self.interpreter.on(verbose=False) 
191   
192          # Replacement user functions. 
193          if user_fns: 
194              for name in user_fns: 
195                  setattr(self.interpreter, name, user_fns[name]) 
196   
197          # Execute the protocol. 
198          try: 
199              # Loop over the models. 
200              for self.diff_model in self.diff_model_list: 
201                  # Wait a little while between diffusion models. 
202                  sleep(1) 
203   
204                  # Set the global model name. 
205                  status.auto_analysis[self.pipe_name].diff_model = self.diff_model 
206   
207                  # Initialise the convergence data structures. 
208                  self.conv_data = Container() 
209                  self.conv_data.chi2 = [] 
210                  self.conv_data.models = [] 
211                  self.conv_data.diff_vals = [] 
212                  if self.diff_model == 'sphere': 
213                      self.conv_data.diff_params = ['tm'] 
214                  elif self.diff_model == 'oblate' or self.diff_model == 'prolate': 
215                      self.conv_data.diff_params = ['tm', 'Da', 'theta', 'phi'] 
216                  elif self.diff_model == 'ellipsoid': 
217                      self.conv_data.diff_params = ['tm', 'Da', 'Dr', 'alpha', 'beta', 'gamma'] 
218                  self.conv_data.spin_ids = [] 
219                  self.conv_data.mf_params = [] 
220                  self.conv_data.mf_vals = [] 
221   
222                  # Execute the analysis for each diffusion model. 
223                  self.execute() 
224   
225          # Clean up. 
226          finally: 
227              # Finish and unlock execution. 
228              status.auto_analysis[self.pipe_name].fin = True 
229              status.current_analysis = None 
230              status.exec_lock.release() 
231   
232   
233 -    def check_vars(self): 
234          """Check that the user has set the variables correctly.""" 
235   
236          # The diff model. 
237          valid_models = ['local_tm', 'sphere', 'oblate', 'prolate', 'ellipsoid', 'final'] 
238          for i in range(len(self.diff_model_list)): 
239              if self.diff_model_list[i] not in valid_models: 
240                  raise RelaxError("The diff_model value '%s' is incorrectly set.  It must be one of %s." % (self.diff_model_list[i], valid_models)) 
241   
242          # Model-free models. 
243          mf_models = ['m0', 'm1', 'm2', 'm3', 'm4', 'm5', 'm6', 'm7', 'm8', 'm9'] 
244          local_tm_models = ['tm0', 'tm1', 'tm2', 'tm3', 'tm4', 'tm5', 'tm6', 'tm7', 'tm8', 'tm9'] 
245          if not isinstance(self.mf_models, list): 
246              raise RelaxError("The self.mf_models user variable must be a list.") 
247          if not isinstance(self.local_tm_models, list): 
248              raise RelaxError("The self.local_tm_models user variable must be a list.") 
249          for i in range(len(self.mf_models)): 
250              if self.mf_models[i] not in mf_models: 
251                  raise RelaxError("The self.mf_models user variable '%s' is incorrectly set.  It must be one of %s." % (self.mf_models, mf_models)) 
252          for i in range(len(self.local_tm_models)): 
253              if self.local_tm_models[i] not in local_tm_models: 
254                  raise RelaxError("The self.local_tm_models user variable '%s' is incorrectly set.  It must be one of %s." % (self.local_tm_models, local_tm_models)) 
255   
256          # Sequence data. 
257          if not exists_mol_res_spin_data(): 
258              raise RelaxNoSequenceError(self.pipe_name) 
259   
260          # Relaxation data. 
261          if not hasattr(cdp, 'ri_ids') or len(cdp.ri_ids) == 0: 
262              raise RelaxNoRiError(ri_id) 
263   
264          # Insufficient data. 
265          if len(cdp.ri_ids) <= 3: 
266              raise RelaxError("Insufficient relaxation data, 4 or more data sets are essential for the execution of this script.") 
267   
268          # Spin vars. 
269          for spin, spin_id in spin_loop(return_id=True): 
270              # Skip deselected spins. 
271              if not spin.select: 
272                  continue 
273   
274              # Print. 
275              print("Checking spin '%s'." % spin_id) 
276   
277              # Test if the bond length has been set. 
278              if not hasattr(spin, 'r') or spin.r == None: 
279                  raise RelaxNoValueError("bond length") 
280   
281              # Test if the CSA value has been set. 
282              if not hasattr(spin, 'csa') or spin.csa == None: 
283                  raise RelaxNoValueError("CSA") 
284   
285              # Test if the heteronucleus type has been set. 
286              if not hasattr(spin, 'heteronuc_type') or spin.heteronuc_type == None: 
287                  raise RelaxNoValueError("heteronucleus type") 
288   
289              # Test if the proton type has been set. 
290              if not hasattr(spin, 'proton_type') or spin.proton_type == None: 
291                  raise RelaxNoValueError("proton type") 
292   
293          # Min vars. 
294          if not isinstance(self.grid_inc, int): 
295              raise RelaxError("The grid_inc user variable '%s' is incorrectly set.  It should be an integer." % self.grid_inc) 
296          if not isinstance(self.diff_tensor_grid_inc, dict): 
297              raise RelaxError("The diff_tensor_grid_inc user variable '%s' is incorrectly set.  It should be a dictionary." % self.diff_tensor_grid_inc) 
298          for tensor in ['sphere', 'prolate', 'oblate', 'ellipsoid']: 
299              if not self.diff_tensor_grid_inc.has_key(tensor): 
300                  raise RelaxError("The diff_tensor_grid_inc user variable '%s' is incorrectly set.  It should contain the '%s' key." % (self.diff_tensor_grid_inc, tensor)) 
301              if not isinstance(self.diff_tensor_grid_inc[tensor], int): 
302                  raise RelaxError("The diff_tensor_grid_inc user variable '%s' is incorrectly set.  The value corresponding to the key '%s' should be an integer." % (self.diff_tensor_grid_inc, tensor)) 
303          if not isinstance(self.min_algor, str): 
304              raise RelaxError("The min_algor user variable '%s' is incorrectly set.  It should be a string." % self.min_algor) 
305          if not isinstance(self.mc_sim_num, int): 
306              raise RelaxError("The mc_sim_num user variable '%s' is incorrectly set.  It should be an integer." % self.mc_sim_num) 
307   
308          # Looping. 
309          if not isinstance(self.conv_loop, bool): 
310              raise RelaxError("The conv_loop user variable '%s' is incorrectly set.  It should be one of the booleans True or False." % self.conv_loop) 
311   
312   
313 -    def convergence(self): 
314          """Test for the convergence of the global model.""" 
315   
316          # Print out. 
317          print("\n\n\n") 
318          print("#####################") 
319          print("# Convergence tests #") 
320          print("#####################\n") 
321   
322          # Maximum number of iterations reached. 
323          if self.max_iter and self.round > self.max_iter: 
324              print("Maximum number of global iterations reached.  Terminating the protocol before convergence has been reached.") 
325              return True 
326   
327          # Store the data of the current data pipe. 
328          self.conv_data.chi2.append(cdp.chi2) 
329   
330          # Create a string representation of the model-free models of the current data pipe. 
331          curr_models = '' 
332          for spin in spin_loop(): 
333              if hasattr(spin, 'model'): 
334                  if not spin.model == 'None': 
335                      curr_models = curr_models + spin.model 
336          self.conv_data.models.append(curr_models) 
337   
338          # Store the diffusion tensor parameters. 
339          self.conv_data.diff_vals.append([]) 
340          for param in self.conv_data.diff_params: 
341              # Get the parameter values. 
342              self.conv_data.diff_vals[-1].append(getattr(cdp.diff_tensor, param)) 
343   
344          # Store the model-free parameters. 
345          self.conv_data.mf_vals.append([]) 
346          self.conv_data.mf_params.append([]) 
347          self.conv_data.spin_ids.append([]) 
348          for spin, spin_id in spin_loop(return_id=True): 
349              # Skip spin systems with no 'params' object. 
350              if not hasattr(spin, 'params'): 
351                  continue 
352   
353              # Add the spin ID, parameters, and empty value list. 
354              self.conv_data.spin_ids[-1].append(spin_id) 
355              self.conv_data.mf_params[-1].append([]) 
356              self.conv_data.mf_vals[-1].append([]) 
357   
358              # Loop over the parameters. 
359              for j in xrange(len(spin.params)): 
360                  # Get the parameters and values. 
361                  self.conv_data.mf_params[-1][-1].append(spin.params[j]) 
362                  self.conv_data.mf_vals[-1][-1].append(getattr(spin, lower(spin.params[j]))) 
363   
364          # No need for tests. 
365          if self.round == 1: 
366              print("First round of optimisation, skipping the convergence tests.\n\n\n") 
367              return False 
368   
369          # Loop over the iterations. 
370          converged = False 
371          for i in range(self.start_round, self.round - 1): 
372              # Print out. 
373              print("\n\n\n# Comparing the current iteration to iteration %i.\n" % (i+1)) 
374   
375              # Index. 
376              index = i - self.start_round 
377   
378              # Chi-squared test. 
379              print("Chi-squared test:") 
380              print("    chi2 (iter %i):  %s" % (i+1, self.conv_data.chi2[index])) 
381              print("        (as an IEEE-754 byte array:  %s)" % floatAsByteArray(self.conv_data.chi2[index])) 
382              print("    chi2 (iter %i):  %s" % (self.round, self.conv_data.chi2[-1])) 
383              print("        (as an IEEE-754 byte array:  %s)" % floatAsByteArray(self.conv_data.chi2[-1])) 
384              print("    chi2 (difference):  %s" % (self.conv_data.chi2[index] - self.conv_data.chi2[-1])) 
385              if self.conv_data.chi2[index] == self.conv_data.chi2[-1]: 
386                  print("    The chi-squared value has converged.\n") 
387              else: 
388                  print("    The chi-squared value has not converged.\n") 
389                  continue 
390   
391              # Identical model-free model test. 
392              print("Identical model-free models test:") 
393              if self.conv_data.models[index] == self.conv_data.models[-1]: 
394                  print("    The model-free models have converged.\n") 
395              else: 
396                  print("    The model-free models have not converged.\n") 
397                  continue 
398   
399              # Identical diffusion tensor parameter value test. 
400              print("Identical diffusion tensor parameter test:") 
401              params_converged = True 
402              for k in range(len(self.conv_data.diff_params)): 
403                  # Test if not identical. 
404                  if self.conv_data.diff_vals[index][k] != self.conv_data.diff_vals[-1][k]: 
405                      print("    Parameter:   %s" % param) 
406                      print("    Value (iter %i):  %s" % (i+1, self.conv_data.diff_vals[index][k])) 
407                      print("        (as an IEEE-754 byte array:  %s)" % floatAsByteArray(self.conv_data.diff_vals[index][k])) 
408                      print("    Value (iter %i):  %s" % (self.round, self.conv_data.diff_vals[-1][k])) 
409                      print("        (as an IEEE-754 byte array:  %s)" % floatAsByteArray(self.conv_data.diff_vals[-1][k])) 
410                      print("    The diffusion parameters have not converged.\n") 
411                      params_converged = False 
412                      break 
413              if not params_converged: 
414                  continue 
415              print("    The diffusion tensor parameters have converged.\n") 
416   
417              # Identical model-free parameter value test. 
418              print("\nIdentical model-free parameter test:") 
419              if len(self.conv_data.spin_ids[index]) != len(self.conv_data.spin_ids[-1]): 
420                  print("    Different number of spins.") 
421                  continue 
422              for j in range(len(self.conv_data.spin_ids[-1])): 
423                  # Loop over the parameters. 
424                  for k in range(len(self.conv_data.mf_params[-1][j])): 
425                      # Test if not identical. 
426                      if self.conv_data.mf_vals[index][j][k] != self.conv_data.mf_vals[-1][j][k]: 
427                          print("    Spin ID:     %s" % self.conv_data.spin_ids[-1][j]) 
428                          print("    Parameter:   %s" % self.conv_data.mf_params[-1][j][k]) 
429                          print("    Value (iter %i): %s" % (i+1, self.conv_data.mf_vals[index][j][k])) 
430                          print("        (as an IEEE-754 byte array:  %s)" % floatAsByteArray(self.conv_data.mf_vals[index][j][k])) 
431                          print("    Value (iter %i): %s" % (self.round, self.conv_data.mf_vals[-1][j][k])) 
432                          print("        (as an IEEE-754 byte array:  %s)" % floatAsByteArray(self.conv_data.mf_vals[index][j][k])) 
433                          print("    The model-free parameters have not converged.\n") 
434                          params_converged = False 
435                          break 
436              if not params_converged: 
437                  continue 
438              print("    The model-free parameters have converged.\n") 
439   
440              # Convergence. 
441              converged = True 
442              break 
443   
444   
445          # Final print out. 
446          ################## 
447   
448          print("\nConvergence:") 
449          if converged: 
450              # Update the status. 
451              status.auto_analysis[self.pipe_name].convergence = True 
452   
453              # Print out. 
454              print("    [ Yes ]") 
455   
456              # Return the termination condition. 
457              return True 
458          else: 
459              # Print out. 
460              print("    [ No ]") 
461   
462              # Return False to not terminate. 
463              return False 
464   
465   
466 -    def determine_rnd(self, model=None): 
467          """Function for returning the name of next round of optimisation.""" 
468   
469          # Get a list of all files in the directory model.  If no directory exists, set the round to 'init' or 0. 
470          try: 
471              dir_list = listdir(self.results_dir+sep+model) 
472          except: 
473              return 0 
474   
475          # Set the round to 'init' or 0 if there is no directory called 'init'. 
476          if 'init' not in dir_list: 
477              return 0 
478   
479          # Create a list of all files which begin with 'round_'. 
480          rnd_dirs = [] 
481          for file in dir_list: 
482              if search('^round_', file): 
483                  rnd_dirs.append(file) 
484   
485          # Create a sorted list of integer round numbers. 
486          numbers = [] 
487          for dir in rnd_dirs: 
488              try: 
489                  numbers.append(int(dir[6:])) 
490              except: 
491                  pass 
492          numbers.sort() 
493   
494          # No directories beginning with 'round_' exist, set the round to 1. 
495          if not len(numbers): 
496              return 1 
497   
498          # The highest number. 
499          max_round = numbers[-1] 
500   
501          # Check that the opt/results file exists for the round (working backwards). 
502          for i in range(max_round, -1, -1): 
503              # Assume the round is complete. 
504              complete_round = i 
505   
506              # The file root. 
507              file_root = self.results_dir + sep + model + sep + "round_%i" % i + sep + 'opt' + sep + 'results' 
508   
509              # Stop looping when the opt/results file is found. 
510              if access(file_root + '.bz2', F_OK): 
511                  break 
512              if access(file_root + '.gz', F_OK): 
513                  break 
514              if access(file_root, F_OK): 
515                  break 
516   
517          # No round, so assume the initial state. 
518          if complete_round == 0: 
519              return 0 
520   
521          # Determine the number for the next round (add 1 to the highest completed round). 
522          return complete_round + 1 
523   
524   
525 -    def execute(self): 
526          """Execute the protocol.""" 
527   
528          # MI - Local tm. 
529          ################ 
530   
531          if self.diff_model == 'local_tm': 
532              # Base directory to place files into. 
533              self.base_dir = self.results_dir+'local_tm'+sep 
534   
535              # Sequential optimisation of all model-free models (function must be modified to suit). 
536              self.multi_model(local_tm=True) 
537   
538              # Model selection. 
539              self.model_selection(modsel_pipe='aic', dir=self.base_dir + 'aic') 
540   
541   
542          # Diffusion models MII to MV. 
543          ############################# 
544   
545          elif self.diff_model == 'sphere' or self.diff_model == 'prolate' or self.diff_model == 'oblate' or self.diff_model == 'ellipsoid': 
546              # No local_tm directory! 
547              dir_list = listdir(self.results_dir) 
548              if 'local_tm' not in dir_list: 
549                  raise RelaxError("The local_tm model must be optimised first.") 
550   
551              # The initial round of optimisation - not zero if calculations were interrupted. 
552              self.start_round = self.determine_rnd(model=self.diff_model) 
553   
554              # Loop until convergence if conv_loop is set, otherwise just loop once. 
555              # This looping could be made much cleaner by removing the dependence on the determine_rnd() function. 
556              while True: 
557                  # Determine which round of optimisation to do (init, round_1, round_2, etc). 
558                  self.round = self.determine_rnd(model=self.diff_model) 
559                  status.auto_analysis[self.pipe_name].round = self.round 
560   
561                  # Inital round of optimisation for diffusion models MII to MV. 
562                  if self.round == 0: 
563                      # Base directory to place files into. 
564                      self.base_dir = self.results_dir+self.diff_model+sep+'init'+sep 
565   
566                      # Run name. 
567                      name = self.diff_model 
568   
569                      # Create the data pipe (deleting the old one if it exists). 
570                      if has_pipe(name): 
571                          self.interpreter.pipe.delete(name) 
572                      self.interpreter.pipe.create(name, 'mf') 
573   
574                      # Load the local tm diffusion model MI results. 
575                      self.interpreter.results.read(file='results', dir=self.results_dir+'local_tm'+sep+'aic') 
576   
577                      # Remove the tm parameter. 
578                      self.interpreter.model_free.remove_tm() 
579   
580                      # Add an arbitrary diffusion tensor which will be optimised. 
581                      if self.diff_model == 'sphere': 
582                          self.interpreter.diffusion_tensor.init(10e-9, fixed=False) 
583                          inc = self.diff_tensor_grid_inc['sphere'] 
584                      elif self.diff_model == 'prolate': 
585                          self.interpreter.diffusion_tensor.init((10e-9, 0, 0, 0), spheroid_type='prolate', fixed=False) 
586                          inc = self.diff_tensor_grid_inc['prolate'] 
587                      elif self.diff_model == 'oblate': 
588                          self.interpreter.diffusion_tensor.init((10e-9, 0, 0, 0), spheroid_type='oblate', fixed=False) 
589                          inc = self.diff_tensor_grid_inc['oblate'] 
590                      elif self.diff_model == 'ellipsoid': 
591                          self.interpreter.diffusion_tensor.init((10e-09, 0, 0, 0, 0, 0), fixed=False) 
592                          inc = self.diff_tensor_grid_inc['ellipsoid'] 
593   
594                      # Minimise just the diffusion tensor. 
595                      self.interpreter.fix('all_spins') 
596                      self.interpreter.grid_search(inc=inc) 
597                      self.interpreter.minimise(self.min_algor, func_tol=self.opt_func_tol, max_iterations=self.opt_max_iterations) 
598   
599                      # Write the results. 
600                      self.interpreter.results.write(file='results', dir=self.base_dir, force=True) 
601   
602   
603                  # Normal round of optimisation for diffusion models MII to MV. 
604                  else: 
605                      # Base directory to place files into. 
606                      self.base_dir = self.results_dir+self.diff_model + sep+'round_'+repr(self.round)+sep 
607   
608                      # Load the optimised diffusion tensor from either the previous round. 
609                      self.load_tensor() 
610   
611                      # Sequential optimisation of all model-free models (function must be modified to suit). 
612                      self.multi_model() 
613   
614                      # Model selection. 
615                      self.model_selection(modsel_pipe='aic', dir=self.base_dir + 'aic') 
616   
617                      # Final optimisation of all diffusion and model-free parameters. 
618                      self.interpreter.fix('all', fixed=False) 
619   
620                      # Minimise all parameters. 
621                      self.interpreter.minimise(self.min_algor, func_tol=self.opt_func_tol, max_iterations=self.opt_max_iterations) 
622   
623                      # Write the results. 
624                      dir = self.base_dir + 'opt' 
625                      self.interpreter.results.write(file='results', dir=dir, force=True) 
626   
627                      # Test for convergence. 
628                      converged = self.convergence() 
629   
630                      # Break out of the infinite while loop if automatic looping is not activated or if convergence has occurred. 
631                      if converged or not self.conv_loop: 
632                          break 
633   
634              # Unset the status. 
635              status.auto_analysis[self.pipe_name].round = None 
636   
637   
638          # Final run. 
639          ############ 
640   
641          elif self.diff_model == 'final': 
642              # Diffusion model selection. 
643              ############################ 
644   
645              # All the global diffusion models to be used in the model selection. 
646              self.pipes = ['local_tm', 'sphere', 'prolate', 'oblate', 'ellipsoid'] 
647   
648              # Remove all temporary pipes used in this auto-analysis. 
649              for name in pipe_names(): 
650                  if name in self.pipes + self.mf_models + self.local_tm_models + ['aic', 'previous']: 
651                      self.interpreter.pipe.delete(name) 
652   
653              # Missing optimised model. 
654              dir_list = listdir(self.results_dir) 
655              for name in self.pipes: 
656                  if name not in dir_list: 
657                      raise RelaxError("The %s model must be optimised first." % name) 
658   
659              # Create the local_tm data pipe. 
660              self.interpreter.pipe.create('local_tm', 'mf') 
661   
662              # Load the local tm diffusion model MI results. 
663              self.interpreter.results.read(file='results', dir=self.results_dir+'local_tm'+sep+'aic') 
664   
665              # Loop over models MII to MV. 
666              for model in ['sphere', 'prolate', 'oblate', 'ellipsoid']: 
667                  # Determine which was the last round of optimisation for each of the models. 
668                  self.round = self.determine_rnd(model=model) - 1 
669   
670                  # If no directories begining with 'round_' exist, the script has not been properly utilised! 
671                  if self.round < 1: 
672                      # Construct the name of the diffusion tensor. 
673                      name = model 
674                      if model == 'prolate' or model == 'oblate': 
675                          name = name + ' spheroid' 
676   
677                      # Throw an error to prevent misuse of the script. 
678                      raise RelaxError("Multiple rounds of optimisation of the " + name + " (between 8 to 15) are required for the proper execution of this script.") 
679   
680                  # Create the data pipe. 
681                  self.interpreter.pipe.create(model, 'mf') 
682   
683                  # Load the diffusion model results. 
684                  self.interpreter.results.read(file='results', dir=self.results_dir+model + sep+'round_'+repr(self.round)+sep+'opt') 
685   
686              # Model selection between MI to MV. 
687              self.model_selection(modsel_pipe='final', write_flag=False) 
688   
689   
690              # Monte Carlo simulations. 
691              ########################## 
692   
693              # Fix the diffusion tensor, if it exists. 
694              if hasattr(get_pipe('final'), 'diff_tensor'): 
695                  self.interpreter.fix('diff') 
696   
697              # Simulations. 
698              self.interpreter.monte_carlo.setup(number=self.mc_sim_num) 
699              self.interpreter.monte_carlo.create_data() 
700              self.interpreter.monte_carlo.initial_values() 
701              self.interpreter.minimise(self.min_algor, func_tol=self.opt_func_tol, max_iterations=self.opt_max_iterations) 
702              self.interpreter.eliminate() 
703              self.interpreter.monte_carlo.error_analysis() 
704   
705   
706              # Write the final results. 
707              ########################## 
708   
709              # Create results files and plots of the data. 
710              self.write_results() 
711   
712   
713          # Unknown script behaviour. 
714          ########################### 
715   
716          else: 
717              raise RelaxError("Unknown diffusion model, change the value of 'self.diff_model'") 
718   
719   
720 -    def load_tensor(self): 
721          """Function for loading the optimised diffusion tensor.""" 
722   
723          # Create the data pipe for the previous data (deleting the old data pipe first if necessary). 
724          if has_pipe('previous'): 
725              self.interpreter.pipe.delete('previous') 
726          self.interpreter.pipe.create('previous', 'mf') 
727   
728          # Load the optimised diffusion tensor from the initial round. 
729          if self.round == 1: 
730              self.interpreter.results.read('results', self.results_dir+self.diff_model + sep+'init') 
731   
732          # Load the optimised diffusion tensor from the previous round. 
733          else: 
734              self.interpreter.results.read('results', self.results_dir+self.diff_model + sep+'round_'+repr(self.round-1)+sep+'opt') 
735   
736   
737 -    def model_selection(self, modsel_pipe=None, dir=None, write_flag=True): 
738          """Model selection function.""" 
739   
740          # Model selection (delete the model selection pipe if it already exists). 
741          if has_pipe(modsel_pipe): 
742              self.interpreter.pipe.delete(modsel_pipe) 
743          self.interpreter.model_selection(method='AIC', modsel_pipe=modsel_pipe, pipes=self.pipes) 
744   
745          # Write the results. 
746          if write_flag: 
747              self.interpreter.results.write(file='results', dir=dir, force=True) 
748   
749   
750 -    def multi_model(self, local_tm=False): 
751          """Function for optimisation of all model-free models.""" 
752   
753          # Set the data pipe names (also the names of preset model-free models). 
754          if local_tm: 
755              self.pipes = self.local_tm_models 
756          else: 
757              self.pipes = self.mf_models 
758   
759          # Loop over the data pipes. 
760          for name in self.pipes: 
761              # Place the model name into the status container. 
762              status.auto_analysis[self.pipe_name].current_model = name 
763   
764              # Create the data pipe (by copying). 
765              if has_pipe(name): 
766                  self.interpreter.pipe.delete(name) 
767              self.interpreter.pipe.copy(self.pipe_name, name) 
768              self.interpreter.pipe.switch(name) 
769   
770              # Copy the diffusion tensor from the 'opt' data pipe and prevent it from being minimised. 
771              if not local_tm: 
772                  self.interpreter.diffusion_tensor.copy('previous') 
773                  self.interpreter.fix('diff') 
774   
775              # Select the model-free model. 
776              self.interpreter.model_free.select_model(model=name) 
777   
778              # Minimise. 
779              self.interpreter.grid_search(inc=self.grid_inc) 
780              self.interpreter.minimise(self.min_algor, func_tol=self.opt_func_tol, max_iterations=self.opt_max_iterations) 
781   
782              # Model elimination. 
783              self.interpreter.eliminate() 
784   
785              # Write the results. 
786              dir = self.base_dir + name 
787              self.interpreter.results.write(file='results', dir=dir, force=True) 
788   
789          # Unset the status. 
790          status.auto_analysis[self.pipe_name].current_model = None 
791   
792   
793 -    def status_setup(self): 
794          """Initialise the status object.""" 
795   
796          # Initialise the status object for this auto-analysis. 
797          status.init_auto_analysis(self.pipe_name, type='dauvergne_protocol') 
798          status.current_analysis = self.pipe_name 
799   
800          # The global diffusion model. 
801          status.auto_analysis[self.pipe_name].diff_model = None 
802   
803          # The round of optimisation, i.e. the global iteration. 
804          status.auto_analysis[self.pipe_name].round = None 
805   
806          # The list of model-free local tm models for optimisation, i.e. the global iteration. 
807          status.auto_analysis[self.pipe_name].local_tm_models = self.local_tm_models 
808   
809          # The list of model-free models for optimisation, i.e. the global iteration. 
810          status.auto_analysis[self.pipe_name].mf_models = self.mf_models 
811   
812          # The current model-free model. 
813          status.auto_analysis[self.pipe_name].current_model = None 
814   
815          # The maximum number of iterations of the global model. 
816          status.auto_analysis[self.pipe_name].max_iter = self.max_iter 
817   
818          # The convergence of the global model. 
819          status.auto_analysis[self.pipe_name].convergence = False 
820   
821   
822 -    def write_results(self): 
823          """Create Grace plots of the final model-free results.""" 
824   
825          # Save the results file. 
826          dir = self.results_dir + 'final' 
827          self.interpreter.results.write(file='results', dir=dir, force=True) 
828   
829          # The Grace plots. 
830          dir = self.results_dir + 'final' + sep + 'grace' 
831          self.interpreter.grace.write(x_data_type='spin', y_data_type='s2',  file='s2.agr',        dir=dir, force=True) 
832          self.interpreter.grace.write(x_data_type='spin', y_data_type='s2f', file='s2f.agr',       dir=dir, force=True) 
833          self.interpreter.grace.write(x_data_type='spin', y_data_type='s2s', file='s2s.agr',       dir=dir, force=True) 
834          self.interpreter.grace.write(x_data_type='spin', y_data_type='te',  file='te.agr',        dir=dir, force=True) 
835          self.interpreter.grace.write(x_data_type='spin', y_data_type='tf',  file='tf.agr',        dir=dir, force=True) 
836          self.interpreter.grace.write(x_data_type='spin', y_data_type='ts',  file='ts.agr',        dir=dir, force=True) 
837          self.interpreter.grace.write(x_data_type='spin', y_data_type='rex', file='rex.agr',       dir=dir, force=True) 
838          self.interpreter.grace.write(x_data_type='s2',   y_data_type='te',  file='s2_vs_te.agr',  dir=dir, force=True) 
839          self.interpreter.grace.write(x_data_type='s2',   y_data_type='rex', file='s2_vs_rex.agr', dir=dir, force=True) 
840          self.interpreter.grace.write(x_data_type='te',   y_data_type='rex', file='te_vs_rex.agr', dir=dir, force=True) 
841   
842          # Write the values to text files. 
843          dir = self.results_dir + 'final' 
844          self.interpreter.value.write(param='s2',       file='s2.txt',       dir=dir, force=True) 
845          self.interpreter.value.write(param='s2f',      file='s2f.txt',      dir=dir, force=True) 
846          self.interpreter.value.write(param='s2s',      file='s2s.txt',      dir=dir, force=True) 
847          self.interpreter.value.write(param='te',       file='te.txt',       dir=dir, force=True) 
848          self.interpreter.value.write(param='tf',       file='tf.txt',       dir=dir, force=True) 
849          self.interpreter.value.write(param='ts',       file='ts.txt',       dir=dir, force=True) 
850          self.interpreter.value.write(param='rex',      file='rex.txt',      dir=dir, force=True) 
851          self.interpreter.value.write(param='local_tm', file='local_tm.txt', dir=dir, force=True) 
852   
853          # Create the PyMOL macros. 
854          dir = self.results_dir + 'final' + sep + 'pymol' 
855          self.interpreter.pymol.macro_write(data_type='s2',        dir=dir, force=True) 
856          self.interpreter.pymol.macro_write(data_type='s2f',       dir=dir, force=True) 
857          self.interpreter.pymol.macro_write(data_type='s2s',       dir=dir, force=True) 
858          self.interpreter.pymol.macro_write(data_type='amp_fast',  dir=dir, force=True) 
859          self.interpreter.pymol.macro_write(data_type='amp_slow',  dir=dir, force=True) 
860          self.interpreter.pymol.macro_write(data_type='te',        dir=dir, force=True) 
861          self.interpreter.pymol.macro_write(data_type='tf',        dir=dir, force=True) 
862          self.interpreter.pymol.macro_write(data_type='ts',        dir=dir, force=True) 
863          self.interpreter.pymol.macro_write(data_type='time_fast', dir=dir, force=True) 
864          self.interpreter.pymol.macro_write(data_type='time_slow', dir=dir, force=True) 
865          self.interpreter.pymol.macro_write(data_type='rex',       dir=dir, force=True) 
866   
867          # Create the Molmol macros. 
868          dir = self.results_dir + 'final' + sep + 'molmol' 
869          self.interpreter.molmol.macro_write(data_type='s2',        dir=dir, force=True) 
870          self.interpreter.molmol.macro_write(data_type='s2f',       dir=dir, force=True) 
871          self.interpreter.molmol.macro_write(data_type='s2s',       dir=dir, force=True) 
872          self.interpreter.molmol.macro_write(data_type='amp_fast',  dir=dir, force=True) 
873          self.interpreter.molmol.macro_write(data_type='amp_slow',  dir=dir, force=True) 
874          self.interpreter.molmol.macro_write(data_type='te',        dir=dir, force=True) 
875          self.interpreter.molmol.macro_write(data_type='tf',        dir=dir, force=True) 
876          self.interpreter.molmol.macro_write(data_type='ts',        dir=dir, force=True) 
877          self.interpreter.molmol.macro_write(data_type='time_fast', dir=dir, force=True) 
878          self.interpreter.molmol.macro_write(data_type='time_slow', dir=dir, force=True) 
879          self.interpreter.molmol.macro_write(data_type='rex',       dir=dir, force=True) 
880   
881          # Create a diffusion tensor representation of the tensor. 
882          dir = self.results_dir + 'final' 
883          self.interpreter.structure.create_diff_tensor_pdb(file="tensor.pdb", dir=dir, force=True) 
884   
885   
886   
887 -class Container: 
888      """Empty container for data storage.""" 
889