1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 """The automatic relaxation curve fitting protocol."""
25
26
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
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
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
70 status.exec_lock.acquire(pipe_bundle, mode='auto-analysis')
71
72
73 status.init_auto_analysis(pipe_bundle, type='relax_fit')
74 status.current_analysis = pipe_bundle
75
76
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
90 self.check_vars()
91
92
93 if self.pipe_name != cdp_name():
94 switch(self.pipe_name)
95
96
97 self.interpreter = Interpreter(show_script=False, raise_relax_error=True)
98 self.interpreter.populate_self()
99 self.interpreter.on(verbose=False)
100
101
102 self.run()
103
104
105 status.auto_analysis[self.pipe_bundle].fin = True
106 status.current_analysis = None
107 status.exec_lock.release()
108
109
111 """Set up and run the curve-fitting."""
112
113
114 self.error_analysis()
115
116
117 self.interpreter.minimise.grid_search(inc=self.grid_inc)
118
119
120 self.interpreter.minimise.execute('newton', scaling=False, constraints=False)
121
122
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
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
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
145 self.interpreter.results.write(file='results', dir=self.results_dir, force=True)
146
147
148 self.interpreter.grace.write(y_data_type='chi2', file='chi2.agr', dir=self.grace_dir, force=True)
149 self.interpreter.grace.write(y_data_type='i0', file='i0.agr', dir=self.grace_dir, force=True)
150 if iinf:
151 self.interpreter.grace.write(y_data_type='iinf', file='iinf.agr', dir=self.grace_dir, force=True)
152 self.interpreter.grace.write(y_data_type='rx', file=self.file_root+'.agr', dir=self.grace_dir, force=True)
153 self.interpreter.grace.write(x_data_type='relax_times', y_data_type='peak_intensity', file='intensities.agr', dir=self.grace_dir, force=True)
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)
155
156
157
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
163 script_grace2images(file=file)
164
165
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
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
185 self.interpreter.state.save(state=self.file_root+'.save', dir=self.results_dir, force=True)
186
187
189 """Perform an error analysis of the peak intensities."""
190
191
192 section(file=sys.stdout, text="Error analysis", prespace=2)
193
194
195 precalc = True
196 for spin in spin_loop(skip_desel=True):
197
198 if not hasattr(spin, 'peak_intensity_err'):
199 precalc = False
200 break
201
202
203 for id in cdp.spectrum_ids:
204 if id not in spin.peak_intensity_err:
205 precalc = False
206 break
207
208
209 if precalc:
210 print("Skipping the error analysis as it has already been performed.")
211 return
212
213
214
215
216 has_dub = False
217
218 if not hasattr(cdp, 'replicates'):
219
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
227 dublicates = [(val, [i for i in range(len(all_times)) if all_times[i] == val]) for val in all_times]
228
229
230 list_dub_mapping = []
231 for i, dub in enumerate(dublicates):
232
233 cur_spectrum_id = all_id[i]
234
235
236 time, list_index_occur = dub
237
238
239 id_list = []
240 if len(list_index_occur) > 1:
241
242 has_dub = True
243
244 for list_index in list_index_occur:
245 id_list.append(all_id[list_index])
246
247
248 list_dub_mapping.append((cur_spectrum_id, id_list))
249
250
251 if has_dub:
252
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
258 self.interpreter.spectrum.error_analysis()
259
260
262 """Check that the user has set the variables correctly."""
263
264
265 if not has_pipe(self.pipe_name):
266 raise RelaxNoPipeError(self.pipe_name)
267