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