1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25 """The R1 and R2 exponential relaxation curve fitting API object."""
26
27
28 from minfx.generic import generic_minimise
29 from minfx.grid import grid
30 from numpy import asarray, dot, float64, transpose, zeros
31 from numpy.linalg import inv
32 from re import match, search
33 import sys
34 from warnings import warn
35
36
37 from dep_check import C_module_exp_fn
38 from lib.errors import RelaxError, RelaxNoModelError
39 from lib.text.sectioning import subsection
40 from lib.warnings import RelaxDeselectWarning
41 from pipe_control.mol_res_spin import check_mol_res_spin_data, return_spin, spin_loop
42 from specific_analyses.api_base import API_base
43 from specific_analyses.api_common import API_common
44 from specific_analyses.relax_fit.checks import check_model_setup
45 from specific_analyses.relax_fit.optimisation import back_calc
46 from specific_analyses.relax_fit.parameter_object import Relax_fit_params
47 from specific_analyses.relax_fit.parameters import assemble_param_vector, disassemble_param_vector, linear_constraints
48 from target_functions.relax_fit_wrapper import Relax_fit_opt
49
50
52 """Class containing functions for relaxation curve fitting."""
53
54
55 instance = None
56
75
76
78 """Return the Jacobian and weights required for parameter errors via the covariance matrix.
79
80 @keyword model_info: The spin container and the spin ID string from the _model_loop_spin() method.
81 @type model_info: SpinContainer instance, str
82 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
83 @type verbosity: int
84 @return: The Jacobian and weight matrices for the given model.
85 @rtype: numpy rank-2 array, numpy rank-2 array
86 """
87
88
89 spin, spin_id = model_info
90
91
92 if not C_module_exp_fn:
93 raise RelaxError("Relaxation curve fitting is not available. Try compiling the C modules on your platform.")
94
95
96 if not (hasattr(spin, 'rx') and hasattr(spin, 'i0')):
97 raise RelaxError("Spin '%s' does not contain optimised 'rx' and 'i0' values. Try execute: minimise.execute(min_algor='Newton', constraints=False)"%(spin_id))
98
99
100 if hasattr(spin, 'g_count'):
101 if getattr(spin, 'g_count') == 0.0:
102 text = "Spin %s contains a gradient count of 0.0. Is the rx parameter optimised? Try execute: minimise.execute(min_algor='Newton', constraints=False)" %(spin_id)
103 warn(RelaxWarning("%s." % text))
104
105
106 if verbosity >= 1:
107
108 top = 2
109 if verbosity >= 2:
110 top += 2
111 subsection(file=sys.stdout, text="Estimating rx error for spin: %s"%spin_id, prespace=top)
112
113
114 values = []
115 errors = []
116 times = []
117 for key in spin.peak_intensity:
118 values.append(spin.peak_intensity[key])
119 errors.append(spin.peak_intensity_err[key])
120 times.append(cdp.relax_times[key])
121
122
123 values = asarray(values)
124 errors = asarray(errors)
125 times = asarray(times)
126
127
128 param_vector = assemble_param_vector(spin=spin)
129 scaling_list = []
130 for i in range(len(spin.params)):
131 scaling_list.append(1.0)
132
133
134 model = Relax_fit_opt(model=spin.model, num_params=len(param_vector), values=values, errors=errors, relax_times=times, scaling_matrix=scaling_list)
135
136
137 jacobian_matrix_exp = transpose(asarray( model.jacobian(param_vector) ) )
138 weights = 1. / errors**2
139
140
141 return jacobian_matrix_exp, weights
142
143
145 """Create the Monte Carlo peak intensity data.
146
147 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method.
148 @type data_id: str
149 @return: The Monte Carlo simulation data.
150 @rtype: list of floats
151 """
152
153
154 mc_data = {}
155
156
157 spin = return_spin(spin_id=data_id)
158
159
160 if not spin.select:
161 return
162
163
164 if not hasattr(spin, 'peak_intensity'):
165 return
166
167
168 if not hasattr(spin, 'model') or not spin.model:
169 raise RelaxNoModelError
170
171
172 for id in cdp.relax_times:
173
174 value = back_calc(spin=spin, relax_time_id=id)
175
176
177 mc_data[id] = value
178
179
180 return mc_data
181
182
184 """Initialise the spin specific data structures.
185
186 @param data: The spin ID string from the _base_data_loop_spin() method.
187 @type data: str
188 @keyword sim: The Monte Carlo simulation flag, which if true will initialise the simulation data structure.
189 @type sim: bool
190 """
191
192
193 spin = return_spin(spin_id=data)
194
195
196 for name in self.data_names(set='params'):
197
198 list_data = [ 'params' ]
199 if name in list_data:
200 init_data = []
201
202
203 else:
204 init_data = None
205
206
207 if not hasattr(spin, name):
208 setattr(spin, name, init_data)
209
210
212 """Return a vector of parameter names.
213
214 @keyword model_info: The spin container and the spin ID string from the _model_loop_spin() method.
215 @type model_info: SpinContainer instance, str
216 @return: The vector of parameter names.
217 @rtype: list of str
218 """
219
220
221 check_model_setup()
222
223
224 spin, spin_id = model_info
225
226
227 return spin.params
228
229
231 """Return a vector of parameter values.
232
233 @keyword model_info: The spin container and the spin ID string from the _model_loop_spin() method.
234 @type model_info: SpinContainer instance, str
235 @keyword sim_index: The optional Monte Carlo simulation index.
236 @type sim_index: int
237 @return: The vector of parameter values.
238 @rtype: list of str
239 """
240
241
242 spin, spin_id = model_info
243
244
245 return assemble_param_vector(spin=spin, sim_index=sim_index)
246
247
248 - def grid_search(self, lower=None, upper=None, inc=None, scaling_matrix=None, constraints=True, verbosity=1, sim_index=None):
249 """The exponential curve fitting grid search method.
250
251 @keyword lower: The per-model lower bounds of the grid search which must be equal to the number of parameters in the model.
252 @type lower: list of lists of numbers
253 @keyword upper: The per-model upper bounds of the grid search which must be equal to the number of parameters in the model.
254 @type upper: list of lists of numbers
255 @keyword inc: The per-model 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.
256 @type inc: list of lists of int
257 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices.
258 @type scaling_matrix: list of numpy rank-2, float64 array or list of None
259 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used.
260 @type constraints: bool
261 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity.
262 @type verbosity: int
263 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised.
264 @type sim_index: int
265 """
266
267
268 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
269
270
271 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling_matrix=None, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):
272 """Relaxation curve fitting minimisation method.
273
274 @keyword min_algor: The minimisation algorithm to use.
275 @type min_algor: str
276 @keyword min_options: An array of options to be used by the minimisation algorithm.
277 @type min_options: array of str
278 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
279 @type func_tol: None or float
280 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
281 @type grad_tol: None or float
282 @keyword max_iterations: The maximum number of iterations for the algorithm.
283 @type max_iterations: int
284 @keyword constraints: If True, constraints are used during optimisation.
285 @type constraints: bool
286 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices.
287 @type scaling_matrix: list of numpy rank-2, float64 array or list of None
288 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
289 @type verbosity: int
290 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
291 @type sim_index: None or int
292 @keyword lower: The per-model 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.
293 @type lower: list of lists of numbers
294 @keyword upper: The per-model 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.
295 @type upper: list of lists of numbers
296 @keyword inc: The per-model 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.
297 @type inc: list of lists of int
298 """
299
300
301 check_mol_res_spin_data()
302
303
304 model_index = 0
305 for spin, spin_id in self.model_loop():
306
307 if not spin.select:
308 continue
309
310
311 if not hasattr(spin, 'peak_intensity'):
312 continue
313
314
315 param_vector = assemble_param_vector(spin=spin)
316
317
318 if scaling_matrix[model_index] is not None:
319 param_vector = dot(inv(scaling_matrix[model_index]), param_vector)
320
321
322 if constraints:
323 A, b = linear_constraints(spin=spin, scaling_matrix=scaling_matrix[model_index])
324 else:
325 A, b = None, None
326
327
328 if verbosity >= 1:
329
330 if verbosity >= 2:
331 print("\n\n")
332
333 string = "Fitting to spin " + repr(spin_id)
334 print("\n\n" + string)
335 print(len(string) * '~')
336
337
338
339
340
341
342 values = []
343 errors = []
344 times = []
345 for key in spin.peak_intensity:
346
347 if sim_index == None:
348 values.append(spin.peak_intensity[key])
349 else:
350 values.append(spin.peak_intensity_sim[sim_index][key])
351
352
353 errors.append(spin.peak_intensity_err[key])
354
355
356 times.append(cdp.relax_times[key])
357
358
359 scaling_list = []
360 if scaling_matrix[model_index] is None:
361 for i in range(len(param_vector)):
362 scaling_list.append(1.0)
363 else:
364 for i in range(len(scaling_matrix[model_index])):
365 scaling_list.append(scaling_matrix[model_index][i, i])
366
367
368 model = Relax_fit_opt(model=spin.model, num_params=len(spin.params), values=values, errors=errors, relax_times=times, scaling_matrix=scaling_list)
369
370
371
372
373
374 if constraints and not match('^[Gg]rid', min_algor):
375 algor = min_options[0]
376 else:
377 algor = min_algor
378
379
380
381
382
383 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor):
384
385 lm_error = zeros(len(spin.relax_times), float64)
386 index = 0
387 for k in range(len(spin.relax_times)):
388 lm_error[index:index+len(relax_error[k])] = relax_error[k]
389 index = index + len(relax_error[k])
390
391 min_options = min_options + (self.relax_fit.lm_dri, lm_error)
392
393
394
395
396
397
398 if search('^[Gg]rid', min_algor):
399 results = grid(func=model.func, args=(), num_incs=inc[model_index], lower=lower[model_index], upper=upper[model_index], A=A, b=b, verbosity=verbosity)
400
401
402 param_vector, chi2, iter_count, warning = results
403 f_count = iter_count
404 g_count = 0.0
405 h_count = 0.0
406
407
408 else:
409 results = generic_minimise(func=model.func, dfunc=model.dfunc, d2func=model.d2func, 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)
410
411
412 if results == None:
413 return
414 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results
415
416
417 if scaling_matrix[model_index] is not None:
418 param_vector = dot(scaling_matrix[model_index], param_vector)
419
420
421 disassemble_param_vector(param_vector=param_vector, spin=spin, sim_index=sim_index)
422
423
424 if sim_index != None:
425
426 spin.chi2_sim[sim_index] = chi2
427
428
429 spin.iter_sim[sim_index] = iter_count
430
431
432 spin.f_count_sim[sim_index] = f_count
433
434
435 spin.g_count_sim[sim_index] = g_count
436
437
438 spin.h_count_sim[sim_index] = h_count
439
440
441 spin.warning_sim[sim_index] = warning
442
443
444
445 else:
446
447 spin.chi2 = chi2
448
449
450 spin.iter = iter_count
451
452
453 spin.f_count = f_count
454
455
456 spin.g_count = g_count
457
458
459 spin.h_count = h_count
460
461
462 spin.warning = warning
463
464
465 model_index += 1
466
467
469 """Deselect spins which have insufficient data to support minimisation.
470
471 @keyword data_check: A flag to signal if the presence of base data is to be checked for.
472 @type data_check: bool
473 @keyword verbose: A flag which if True will allow printouts.
474 @type verbose: bool
475 """
476
477
478 if verbose:
479 print("\nOver-fit spin deselection:")
480
481
482 check_mol_res_spin_data()
483 check_model_setup()
484
485
486 deselect_flag = False
487 for spin, spin_id in spin_loop(return_id=True):
488
489 if not spin.select:
490 continue
491
492
493 if not hasattr(spin, 'peak_intensity'):
494 warn(RelaxDeselectWarning(spin_id, 'missing intensity data'))
495 spin.select = False
496 deselect_flag = True
497 continue
498
499
500 elif len(spin.peak_intensity) < 3:
501 warn(RelaxDeselectWarning(spin_id, 'insufficient data, 3 or more data points are required'))
502 spin.select = False
503 deselect_flag = True
504 continue
505
506
507 if len(spin.peak_intensity) != len(cdp.relax_times):
508 raise RelaxError("The %s peak intensity points of the spin '%s' does not match the expected number of %s (the IDs %s do not match %s)." % (len(spin.peak_intensity), spin_id, len(cdp.relax_times), sorted(spin.peak_intensity.keys()), sorted(cdp.relax_times.keys())))
509
510
511 if verbose and not deselect_flag:
512 print("No spins have been deselected.")
513
514
516 """Function for returning the peak intensity data structure.
517
518 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method.
519 @type data_id: str
520 @return: The peak intensity data structure.
521 @rtype: list of float
522 """
523
524
525 spin = return_spin(spin_id=data_id)
526
527
528 return spin.peak_intensity
529
530
532 """Return the standard deviation data structure.
533
534 @param data_id: The spin identification string, as yielded by the base_data_loop() generator
535 method.
536 @type data_id: str
537 @return: The standard deviation data structure.
538 @rtype: list of float
539 """
540
541
542 spin = return_spin(spin_id=data_id)
543
544
545 return spin.peak_intensity_err
546
547
549 """Pack the Monte Carlo simulation data.
550
551 @param data_id: The spin identification string, as yielded by the base_data_loop() generator method.
552 @type data_id: str
553 @param sim_data: The Monte Carlo simulation data.
554 @type sim_data: list of float
555 """
556
557
558 spin = return_spin(spin_id=data_id)
559
560
561 spin.peak_intensity_sim = sim_data
562