Source Code for Module auto_analyses.relax_fit

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2004-2015 Edward d'Auvergne                                   # 
  4  # Copyright (C) 2014 Troels E. Linnet                                         # 
  5  #                                                                             # 
  6  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  7  #                                                                             # 
  8  # This program is free software: you can redistribute it and/or modify        # 
  9  # it under the terms of the GNU General Public License as published by        # 
 10  # the Free Software Foundation, either version 3 of the License, or           # 
 11  # (at your option) any later version.                                         # 
 12  #                                                                             # 
 13  # This program is distributed in the hope that it will be useful,             # 
 14  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 16  # GNU General Public License for more details.                                # 
 17  #                                                                             # 
 18  # You should have received a copy of the GNU General Public License           # 
 19  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23  # Module docstring. 
 24  """The automatic relaxation curve fitting protocol.""" 
 25   
 26  # Python module imports. 
 27  from os import chmod, sep 
 28  from os.path import expanduser 
 29  from stat import S_IRWXU, S_IRGRP, S_IROTH 
 30  import sys 
 31   
 32  # relax module imports. 
 33  from lib.io import get_file_path, open_write_file 
 34  from lib.plotting.grace import script_grace2images 
 35  from lib.text.sectioning import section 
 36  from pipe_control.mol_res_spin import spin_loop 
 37  from pipe_control.pipes import cdp_name, has_pipe, switch 
 38  from prompt.interpreter import Interpreter 
 39  from status import Status; status = Status() 
 40   
 41   
 42 -class Relax_fit: 
 43 -    def __init__(self, pipe_name=None, pipe_bundle=None, file_root='rx', results_dir=None, grid_inc=11, mc_sim_num=500, view_plots=True): 
 44          """Perform relaxation curve fitting. 
 45   
 46          To use this auto-analysis, a data pipe with all the required data needs to be set up.  This data pipe should contain the following: 
 47   
 48              - All the spins loaded. 
 49              - Unresolved spins deselected. 
 50              - All the peak intensities loaded and relaxation delay times set. 
 51              - Either the baseplane noise RMSD values should be set or replicated spectra loaded. 
 52   
 53          @keyword pipe_name:     The name of the data pipe containing all of the data for the analysis. 
 54          @type pipe_name:        str 
 55          @keyword pipe_bundle:   The data pipe bundle to associate all spawned data pipes with. 
 56          @type pipe_bundle:      str 
 57          @keyword file_root:     File root of the output filea. 
 58          @type file_root:        str 
 59          @keyword results_dir:   The directory where results files are saved. 
 60          @type results_dir:      str 
 61          @keyword grid_inc:      Number of grid search increments. 
 62          @type grid_inc:         int 
 63          @keyword mc_sim_num:    The number of Monte Carlo simulations to be used for error analysis at the end of the analysis. 
 64          @type mc_sim_num:       int 
 65          @keyword view_plots:    Flag to automatically view grace plots after calculation. 
 66          @type view_plots:       bool 
 67          """ 
 68   
 69          # Execution lock. 
 70          status.exec_lock.acquire(pipe_bundle, mode='auto-analysis') 
 71   
 72          # Set up the analysis status object. 
 73          status.init_auto_analysis(pipe_bundle, type='relax_fit') 
 74          status.current_analysis = pipe_bundle 
 75   
 76          # Store the args. 
 77          self.pipe_name = pipe_name 
 78          self.pipe_bundle = pipe_bundle 
 79          self.file_root = file_root 
 80          self.results_dir = results_dir 
 81          if self.results_dir: 
 82              self.grace_dir = results_dir + sep + 'grace' 
 83          else: 
 84              self.grace_dir = 'grace' 
 85          self.mc_sim_num = mc_sim_num 
 86          self.grid_inc = grid_inc 
 87          self.view_plots = view_plots 
 88   
 89          # Data checks. 
 90          self.check_vars() 
 91   
 92          # Set the data pipe to the current data pipe. 
 93          if self.pipe_name != cdp_name(): 
 94              switch(self.pipe_name) 
 95   
 96          # Load the interpreter. 
 97          self.interpreter = Interpreter(show_script=False, raise_relax_error=True) 
 98          self.interpreter.populate_self() 
 99          self.interpreter.on(verbose=False) 
100   
101          # Execute. 
102          self.run() 
103   
104          # Finish and unlock execution. 
105          status.auto_analysis[self.pipe_bundle].fin = True 
106          status.current_analysis = None 
107          status.exec_lock.release() 
108   
109   
110 -    def run(self): 
111          """Set up and run the curve-fitting.""" 
112   
113          # Peak intensity error analysis. 
114          self.error_analysis() 
115   
116          # Grid search. 
117          self.interpreter.minimise.grid_search(inc=self.grid_inc) 
118   
119          # Minimise. 
120          self.interpreter.minimise.execute('newton', scaling=False, constraints=False) 
121   
122          # Monte Carlo simulations. 
123          self.interpreter.monte_carlo.setup(number=self.mc_sim_num) 
124          self.interpreter.monte_carlo.create_data() 
125          self.interpreter.monte_carlo.initial_values() 
126          self.interpreter.minimise.execute('newton', scaling=False, constraints=False) 
127          self.interpreter.monte_carlo.error_analysis() 
128   
129          # Determine the normalisation type and if the Iinf parameter exists. 
130          norm_type = 'last' 
131          iinf = True 
132          for spin in spin_loop(skip_desel=True): 
133              if spin.model not in ['sat', 'inv']: 
134                  norm_type = 'first' 
135                  iinf = False 
136                  break 
137   
138          # Save the relaxation rates and other parameter values. 
139          self.interpreter.value.write(param='i0', file='i0.out', dir=self.results_dir, force=True) 
140          if iinf: 
141              self.interpreter.value.write(param='iinf', file='iinf.out', dir=self.results_dir, force=True) 
142          self.interpreter.value.write(param='rx', file=self.file_root+'.out', dir=self.results_dir, force=True) 
143   
144          # Save the results. 
145          self.interpreter.results.write(file='results', dir=self.results_dir, force=True) 
146   
147          # Create Grace plots of the data. 
148          self.interpreter.grace.write(y_data_type='chi2', file='chi2.agr', dir=self.grace_dir, force=True)    # Minimised chi-squared value. 
149          self.interpreter.grace.write(y_data_type='i0', file='i0.agr', dir=self.grace_dir, force=True)    # Initial peak intensity. 
150          if iinf: 
151              self.interpreter.grace.write(y_data_type='iinf', file='iinf.agr', dir=self.grace_dir, force=True)    # Infinite peak intensity. 
152          self.interpreter.grace.write(y_data_type='rx', file=self.file_root+'.agr', dir=self.grace_dir, force=True)    # Relaxation rate. 
153          self.interpreter.grace.write(x_data_type='relax_times', y_data_type='peak_intensity', file='intensities.agr', dir=self.grace_dir, force=True)    # Average peak intensities. 
154          self.interpreter.grace.write(x_data_type='relax_times', y_data_type='peak_intensity', norm_type=norm_type, norm=True, file='intensities_norm.agr', dir=self.grace_dir, force=True)    # Average peak intensities (normalised). 
155   
156          # Write a python "grace to PNG/EPS/SVG..." conversion script. 
157          # Open the file for writing. 
158          file_name = "grace2images.py" 
159          file_path = get_file_path(file_name=file_name, dir=self.grace_dir) 
160          file = open_write_file(file_name=file_name, dir=self.grace_dir, force=True) 
161   
162          # Write the file. 
163          script_grace2images(file=file) 
164   
165          # Close the batch script, then make it executable (expanding any ~ characters). 
166          file.close() 
167   
168          if self.grace_dir: 
169              dir = expanduser(self.grace_dir) 
170              chmod(dir + sep + file_name, S_IRWXU|S_IRGRP|S_IROTH) 
171          else: 
172              file_name = expanduser(file_name) 
173              chmod(file_name, S_IRWXU|S_IRGRP|S_IROTH) 
174   
175   
176          # Display the Grace plots if selected. 
177          if self.view_plots: 
178              self.interpreter.grace.view(file='chi2.agr', dir=self.grace_dir) 
179              self.interpreter.grace.view(file='i0.agr', dir=self.grace_dir) 
180              self.interpreter.grace.view(file=self.file_root+'.agr', dir=self.grace_dir) 
181              self.interpreter.grace.view(file='intensities.agr', dir=self.grace_dir) 
182              self.interpreter.grace.view(file='intensities_norm.agr', dir=self.grace_dir) 
183   
184          # Save the program state. 
185          self.interpreter.state.save(state=self.file_root+'.save', dir=self.results_dir, force=True) 
186   
187   
188 -    def error_analysis(self): 
189          """Perform an error analysis of the peak intensities.""" 
190   
191          # Printout. 
192          section(file=sys.stdout, text="Error analysis", prespace=2) 
193   
194          # Check if intensity errors have already been calculated by the user. 
195          precalc = True 
196          for spin in spin_loop(skip_desel=True): 
197              # No structure. 
198              if not hasattr(spin, 'peak_intensity_err'): 
199                  precalc = False 
200                  break 
201   
202              # Determine if a spectrum ID is missing from the list. 
203              for id in cdp.spectrum_ids: 
204                  if id not in spin.peak_intensity_err: 
205                      precalc = False 
206                      break 
207   
208          # Skip. 
209          if precalc: 
210              print("Skipping the error analysis as it has already been performed.") 
211              return 
212   
213          # Check if there is replicates, and the user has not specified them. 
214   
215          # Set flag for dublicates. 
216          has_dub = False 
217   
218          if not hasattr(cdp, 'replicates'): 
219              # Collect all times, and matching spectrum ID. 
220              all_times = [] 
221              all_id = [] 
222              for spectrum_id in cdp.relax_times: 
223                  all_times.append(cdp.relax_times[spectrum_id]) 
224                  all_id.append(spectrum_id) 
225   
226              # Get the dublicates. 
227              dublicates = [(val, [i for i in range(len(all_times)) if all_times[i] == val]) for val in all_times] 
228   
229              # Loop over the list of the mapping of times and duplications. 
230              list_dub_mapping = [] 
231              for i, dub in enumerate(dublicates): 
232                  # Get current spectum id. 
233                  cur_spectrum_id = all_id[i] 
234   
235                  # Get the tuple of time and indexes of duplications. 
236                  time, list_index_occur = dub 
237   
238                  # Collect mapping of index to id. 
239                  id_list = [] 
240                  if len(list_index_occur) > 1: 
241                      # There exist dublications. 
242                      has_dub = True 
243   
244                      for list_index in list_index_occur: 
245                          id_list.append(all_id[list_index]) 
246   
247                  # Store to list 
248                  list_dub_mapping.append((cur_spectrum_id, id_list)) 
249   
250          # If there is dublication, then assign them. 
251          if has_dub: 
252              # Assign dublicates. 
253              for spectrum_id, dub_pair in list_dub_mapping: 
254                  if len(dub_pair) > 0: 
255                      self.interpreter.spectrum.replicated(spectrum_ids=dub_pair) 
256   
257          # Run the error analysis. 
258          self.interpreter.spectrum.error_analysis() 
259   
260   
261 -    def check_vars(self): 
262          """Check that the user has set the variables correctly.""" 
263   
264          # The pipe name. 
265          if not has_pipe(self.pipe_name): 
266              raise RelaxNoPipeError(self.pipe_name) 
267