1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """The R1 and R2 exponential relaxation curve fitting API object."""
24
25
26 from minfx.generic import generic_minimise
27 from minfx.grid import grid
28 from numpy import dot, float64, zeros
29 from numpy.linalg import inv
30 from re import match, search
31 from warnings import warn
32
33
34 from dep_check import C_module_exp_fn
35 from lib.errors import RelaxError, RelaxNoModelError, RelaxNoSequenceError
36 from lib.warnings import RelaxDeselectWarning
37 from pipe_control.mol_res_spin import exists_mol_res_spin_data, generate_spin_id_unique, return_spin, spin_loop
38 from specific_analyses.api_base import API_base
39 from specific_analyses.api_common import API_common
40 from specific_analyses.relax_fit.optimisation import back_calc, d2func_wrapper, dfunc_wrapper, func_wrapper, grid_search_setup
41 from specific_analyses.relax_fit.parameter_object import Relax_fit_params
42 from specific_analyses.relax_fit.parameters import assemble_param_vector, assemble_scaling_matrix, disassemble_param_vector, linear_constraints
43
44
45 if C_module_exp_fn:
46 from target_functions.relax_fit import setup
47
48
50 """Class containing functions for relaxation curve fitting."""
51
52
53 instance = None
54
72
73
75 """Create the Monte Carlo peak intensity data.
76
77 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method.
78 @type data_id: str
79 @return: The Monte Carlo simulation data.
80 @rtype: list of floats
81 """
82
83
84 mc_data = {}
85
86
87 spin = return_spin(data_id)
88
89
90 if not spin.select:
91 return
92
93
94 if not hasattr(spin, 'peak_intensity'):
95 return
96
97
98 if not hasattr(spin, 'model') or not spin.model:
99 raise RelaxNoModelError
100
101
102 for id in list(cdp.relax_times.keys()):
103
104 value = back_calc(spin=spin, relax_time_id=id)
105
106
107 mc_data[id] = value
108
109
110 return mc_data
111
112
114 """Initialise the spin specific data structures.
115
116 @param data_cont: The spin container.
117 @type data_cont: SpinContainer instance
118 @keyword sim: The Monte Carlo simulation flag, which if true will initialise the simulation data structure.
119 @type sim: bool
120 """
121
122
123 for name in self.data_names(set='params'):
124
125 list_data = [ 'params' ]
126 if name in list_data:
127 init_data = []
128
129
130 else:
131 init_data = None
132
133
134 if not hasattr(data_cont, name):
135 setattr(data_cont, name, init_data)
136
137
138 - def grid_search(self, lower=None, upper=None, inc=None, constraints=True, verbosity=1, sim_index=None):
139 """The exponential curve fitting grid search method.
140
141 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model.
142 @type lower: array of numbers
143 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model.
144 @type upper: array of numbers
145 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model.
146 @type inc: array of int
147 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used.
148 @type constraints: bool
149 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity.
150 @type verbosity: int
151 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised.
152 @type sim_index: int
153 """
154
155
156 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
157
158
159 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling=True, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):
160 """Relaxation curve fitting minimisation method.
161
162 @keyword min_algor: The minimisation algorithm to use.
163 @type min_algor: str
164 @keyword min_options: An array of options to be used by the minimisation algorithm.
165 @type min_options: array of str
166 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
167 @type func_tol: None or float
168 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
169 @type grad_tol: None or float
170 @keyword max_iterations: The maximum number of iterations for the algorithm.
171 @type max_iterations: int
172 @keyword constraints: If True, constraints are used during optimisation.
173 @type constraints: bool
174 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned.
175 @type scaling: bool
176 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
177 @type verbosity: int
178 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
179 @type sim_index: None or int
180 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search.
181 @type lower: array of numbers
182 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search.
183 @type upper: array of numbers
184 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. This argument is only used when doing a grid search.
185 @type inc: array of int
186 """
187
188
189 if not exists_mol_res_spin_data():
190 raise RelaxNoSequenceError
191
192
193 for spin, mol_name, res_num, res_name in spin_loop(full_info=True):
194
195 if not spin.select:
196 continue
197
198
199 if not hasattr(spin, 'peak_intensity'):
200 continue
201
202
203 param_vector = assemble_param_vector(spin=spin)
204
205
206 scaling_matrix = assemble_scaling_matrix(spin=spin, scaling=scaling)
207 if len(scaling_matrix):
208 param_vector = dot(inv(scaling_matrix), param_vector)
209
210
211 if match('^[Gg]rid', min_algor):
212 inc, lower_new, upper_new = grid_search_setup(spin=spin, param_vector=param_vector, lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix)
213
214
215 if constraints:
216 A, b = linear_constraints(spin=spin, scaling_matrix=scaling_matrix)
217 else:
218 A, b = None, None
219
220
221 if verbosity >= 1:
222
223 spin_id = generate_spin_id_unique(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=spin.num, spin_name=spin.name)
224
225
226 if verbosity >= 2:
227 print("\n\n")
228
229 string = "Fitting to spin " + repr(spin_id)
230 print("\n\n" + string)
231 print(len(string) * '~')
232
233
234
235
236
237
238 keys = list(spin.peak_intensity.keys())
239
240
241 values = []
242 errors = []
243 times = []
244 for key in keys:
245
246 if sim_index == None:
247 values.append(spin.peak_intensity[key])
248 else:
249 values.append(spin.peak_intensity_sim[sim_index][key])
250
251
252 errors.append(spin.peak_intensity_err[key])
253
254
255 times.append(cdp.relax_times[key])
256
257
258 scaling_list = []
259 for i in range(len(scaling_matrix)):
260 scaling_list.append(scaling_matrix[i, i])
261
262 setup(num_params=len(spin.params), num_times=len(values), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list)
263
264
265
266
267
268 if constraints and not match('^[Gg]rid', min_algor):
269 algor = min_options[0]
270 else:
271 algor = min_algor
272
273
274
275
276
277 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor):
278
279 lm_error = zeros(len(spin.relax_times), float64)
280 index = 0
281 for k in range(len(spin.relax_times)):
282 lm_error[index:index+len(relax_error[k])] = relax_error[k]
283 index = index + len(relax_error[k])
284
285 min_options = min_options + (self.relax_fit.lm_dri, lm_error)
286
287
288
289
290
291
292 if search('^[Gg]rid', min_algor):
293 results = grid(func=func_wrapper, args=(), num_incs=inc, lower=lower_new, upper=upper_new, A=A, b=b, verbosity=verbosity)
294
295
296 param_vector, chi2, iter_count, warning = results
297 f_count = iter_count
298 g_count = 0.0
299 h_count = 0.0
300
301
302 else:
303 results = generic_minimise(func=func_wrapper, dfunc=dfunc_wrapper, d2func=d2func_wrapper, args=(), x0=param_vector, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, maxiter=max_iterations, A=A, b=b, full_output=True, print_flag=verbosity)
304
305
306 if results == None:
307 return
308 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results
309
310
311 if scaling:
312 param_vector = dot(scaling_matrix, param_vector)
313
314
315 disassemble_param_vector(param_vector=param_vector, spin=spin, sim_index=sim_index)
316
317
318 if sim_index != None:
319
320 spin.chi2_sim[sim_index] = chi2
321
322
323 spin.iter_sim[sim_index] = iter_count
324
325
326 spin.f_count_sim[sim_index] = f_count
327
328
329 spin.g_count_sim[sim_index] = g_count
330
331
332 spin.h_count_sim[sim_index] = h_count
333
334
335 spin.warning_sim[sim_index] = warning
336
337
338
339 else:
340
341 spin.chi2 = chi2
342
343
344 spin.iter = iter_count
345
346
347 spin.f_count = f_count
348
349
350 spin.g_count = g_count
351
352
353 spin.h_count = h_count
354
355
356 spin.warning = warning
357
358
404
405
407 """Function for returning the peak intensity data structure.
408
409 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method.
410 @type data_id: str
411 @return: The peak intensity data structure.
412 @rtype: list of float
413 """
414
415
416 spin = return_spin(data_id)
417
418
419 return spin.peak_intensity
420
421
423 """Return the standard deviation data structure.
424
425 @param data_id: The spin identification string, as yielded by the base_data_loop() generator
426 method.
427 @type data_id: str
428 @return: The standard deviation data structure.
429 @rtype: list of float
430 """
431
432
433 spin = return_spin(data_id)
434
435
436 return spin.peak_intensity_err
437
438
440 """Pack the Monte Carlo simulation data.
441
442 @param data_id: The spin identification string, as yielded by the base_data_loop() generator method.
443 @type data_id: str
444 @param sim_data: The Monte Carlo simulation data.
445 @type sim_data: list of float
446 """
447
448
449 spin = return_spin(data_id)
450
451
452 spin.peak_intensity_sim = sim_data
453