1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """Module for the optimisation of the relaxation dispersion models."""
24
25
26 from minfx.generic import generic_minimise
27 from minfx.grid import grid
28 from numpy import dot, float64, int32, ones, zeros
29 from numpy.linalg import inv
30 from re import match, search
31 import sys
32
33
34 from dep_check import C_module_exp_fn
35 from lib.check_types import is_float
36 from lib.dispersion.two_point import calc_two_point_r2eff, calc_two_point_r2eff_err
37 from lib.errors import RelaxError
38 from lib.text.sectioning import subsection
39 from multi import Memo, Result_command, Slave_command
40 from pipe_control.mol_res_spin import spin_loop
41 from specific_analyses.relax_disp.checks import check_disp_points, check_exp_type, check_exp_type_fixed_time
42 from specific_analyses.relax_disp.data import average_intensity, count_spins, find_intensity_keys, has_exponential_exp_type, has_proton_mmq_cpmg, loop_exp, loop_exp_frq_offset_point, loop_exp_frq_offset_point_time, loop_frq, loop_offset, loop_time, pack_back_calc_r2eff, return_cpmg_frqs, return_offset_data, return_param_key_from_data, return_r1_data, return_r2eff_arrays, return_spin_lock_nu1
43 from specific_analyses.relax_disp.parameters import assemble_param_vector, assemble_scaling_matrix, disassemble_param_vector, linear_constraints, loop_parameters, param_conversion, param_num
44 from specific_analyses.relax_disp.variables import EXP_TYPE_LIST_CPMG, MODEL_CR72, MODEL_CR72_FULL, MODEL_LIST_MMQ, MODEL_LM63, MODEL_M61, MODEL_M61B, MODEL_MP05, MODEL_TAP03, MODEL_TP02
45 from target_functions.relax_disp import Dispersion
46
47
48 if C_module_exp_fn:
49 from target_functions.relax_fit import setup, func, dfunc, d2func, back_calc_I
50
51
53 """Back-calculation of peak intensity for the given relaxation time.
54
55 @keyword spin: The specific spin data container.
56 @type spin: SpinContainer instance
57 @keyword exp_type: The experiment type.
58 @type exp_type: str
59 @keyword frq: The spectrometer frequency.
60 @type frq: float
61 @keyword offset: For R1rho-type data, the spin-lock offset value in ppm.
62 @type offset: None or float
63 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz).
64 @type point: float
65 @return: The back-calculated peak intensities for the given exponential curve.
66 @rtype: numpy rank-1 float array
67 """
68
69
70 if not has_exponential_exp_type():
71 raise RelaxError("Back-calculation is not allowed for the fixed time experiment types.")
72
73
74 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point)
75
76
77 param_vector = assemble_param_vector(spins=[spin], key=param_key)
78
79
80 scaling_matrix = assemble_scaling_matrix(spins=[spin], key=param_key, scaling=False)
81
82
83 values = []
84 errors = []
85 times = []
86 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point):
87
88 values.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time))
89 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True))
90 times.append(time)
91
92
93 scaling_list = []
94 for i in range(len(scaling_matrix)):
95 scaling_list.append(scaling_matrix[i, i])
96
97
98 setup(num_params=len(param_vector), num_times=len(times), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list)
99
100
101 func(param_vector)
102
103
104 results = back_calc_I()
105
106
107 return results
108
109
110 -def back_calc_r2eff(spin=None, spin_id=None, cpmg_frqs=None, spin_lock_nu1=None, store_chi2=False):
111 """Back-calculation of R2eff/R1rho values for the given spin.
112
113 @keyword spin: The specific spin data container.
114 @type spin: SpinContainer instance
115 @keyword spin_id: The spin ID string for the spin container.
116 @type spin_id: str
117 @keyword cpmg_frqs: The CPMG frequencies to use instead of the user loaded values - to enable interpolation.
118 @type cpmg_frqs: list of lists of numpy rank-1 float arrays
119 @keyword spin_lock_nu1: The spin-lock field strengths to use instead of the user loaded values - to enable interpolation.
120 @type spin_lock_nu1: list of lists of numpy rank-1 float arrays
121 @keyword store_chi2: A flag which if True will cause the spin specific chi-squared value to be stored in the spin container.
122 @type store_chi2: bool
123 @return: The back-calculated R2eff/R1rho value for the given spin.
124 @rtype: numpy rank-1 float array
125 """
126
127
128 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H':
129 return
130
131
132 param_vector = assemble_param_vector(spins=[spin])
133
134
135 scaling_matrix = assemble_scaling_matrix(spins=[spin], scaling=False)
136
137
138 fields = [None]
139 field_count = 1
140 if hasattr(cdp, 'spectrometer_frq_count'):
141 fields = cdp.spectrometer_frq_list
142 field_count = cdp.spectrometer_frq_count
143
144
145 values, errors, missing, frqs, frqs_H, exp_types, relax_times = return_r2eff_arrays(spins=[spin], spin_ids=[spin_id], fields=fields, field_count=field_count)
146
147
148 chemical_shifts, offsets, tilt_angles, Delta_omega, w_eff = return_offset_data(spins=[spin], spin_ids=[spin_id], field_count=field_count, fields=spin_lock_nu1)
149 r1 = return_r1_data(spins=[spin], spin_ids=[spin_id], field_count=field_count)
150
151
152 recalc_tau = True
153 if cpmg_frqs == None and spin_lock_nu1 == None:
154 cpmg_frqs = return_cpmg_frqs(ref_flag=False)
155 spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False)
156
157
158 else:
159 recalc_tau = False
160 values = []
161 errors = []
162 missing = []
163 for exp_type, ei in loop_exp(return_indices=True):
164 values.append([])
165 errors.append([])
166 missing.append([])
167 for si in range(1):
168 values[ei].append([])
169 errors[ei].append([])
170 missing[ei].append([])
171 for frq, mi in loop_frq(return_indices=True):
172 values[ei][si].append([])
173 errors[ei][si].append([])
174 missing[ei][si].append([])
175 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True):
176 if exp_type in EXP_TYPE_LIST_CPMG:
177 num = len(cpmg_frqs[ei][mi][oi])
178 else:
179 num = len(spin_lock_nu1[ei][mi][oi])
180 values[ei][si][mi].append(zeros(num, float64))
181 errors[ei][si][mi].append(ones(num, float64))
182 missing[ei][si][mi].append(zeros(num, int32))
183
184
185 model = Dispersion(model=spin.model, num_params=param_num(spins=[spin]), num_spins=1, num_frq=field_count, exp_types=exp_types, values=values, errors=errors, missing=missing, frqs=frqs, frqs_H=frqs_H, cpmg_frqs=cpmg_frqs, spin_lock_nu1=spin_lock_nu1, chemical_shifts=chemical_shifts, offset=offsets, tilt_angles=tilt_angles, r1=r1, relax_times=relax_times, scaling_matrix=scaling_matrix, recalc_tau=recalc_tau)
186
187
188 chi2 = model.func(param_vector)
189
190
191 if store_chi2:
192 spin.chi2 = chi2
193
194
195 return model.back_calc
196
197
199 """Calculate the R2eff values for fixed relaxation time period data."""
200
201
202 check_exp_type()
203 check_disp_points()
204 check_exp_type_fixed_time()
205
206
207 print("Calculating the R2eff/R1rho values for fixed relaxation time period data.")
208
209
210 for spin, spin_id in spin_loop(return_id=True, skip_desel=True):
211
212 print("Spin '%s'." % spin_id)
213
214
215 if not hasattr(spin, 'peak_intensity'):
216 continue
217
218
219 if not hasattr(spin, 'r2eff'):
220 spin.r2eff = {}
221 if not hasattr(spin, 'r2eff_err'):
222 spin.r2eff_err = {}
223
224
225 for exp_type, frq, offset, point, time in loop_exp_frq_offset_point_time():
226
227 ref_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=None, time=time)
228 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)
229 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point)
230
231
232 missing = False
233 for i in range(len(ref_keys)):
234 if ref_keys[i] not in spin.peak_intensity:
235 missing = True
236 for i in range(len(int_keys)):
237 if int_keys[i] not in spin.peak_intensity:
238 missing = True
239 if missing:
240 continue
241
242
243 ref_intensity = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=None, time=time)
244 ref_intensity_err = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=None, time=time, error=True)
245
246
247 intensity = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)
248 intensity_err = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)
249
250
251 spin.r2eff[param_key] = calc_two_point_r2eff(relax_time=time, I_ref=ref_intensity, I=intensity)
252
253
254 spin.r2eff_err[param_key] = calc_two_point_r2eff_err(relax_time=time, I_ref=ref_intensity, I=intensity, I_ref_err=ref_intensity_err, I_err=intensity_err)
255
256
257 -def grid_search_setup(spins=None, spin_ids=None, param_vector=None, lower=None, upper=None, inc=None, scaling_matrix=None):
258 """The grid search setup function.
259
260 @keyword spins: The list of spin data containers for the block.
261 @type spins: list of SpinContainer instances
262 @keyword spin_ids: The corresponding spin ID strings.
263 @type spin_ids: list of str
264 @keyword param_vector: The parameter vector.
265 @type param_vector: numpy array
266 @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.
267 @type lower: array of numbers
268 @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.
269 @type upper: array of numbers
270 @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.
271 @type inc: array of int
272 @keyword scaling_matrix: The scaling matrix.
273 @type scaling_matrix: numpy diagonal matrix
274 @return: A tuple of the grid size and the minimisation options. For the minimisation options, the first dimension corresponds to the model parameter. The second dimension is a list of the number of increments, the lower bound, and upper bound.
275 @rtype: (int, list of lists [int, float, float])
276 """
277
278
279 n = len(param_vector)
280
281
282 if n == 0:
283 raise RelaxError("Cannot run a grid search on a model with zero parameters.")
284
285
286 if lower != None and len(lower) != n:
287 raise RelaxLenError('lower bounds', n)
288
289
290 if upper != None and len(upper) != n:
291 raise RelaxLenError('upper bounds', n)
292
293
294 if isinstance(inc, list) and len(inc) != n:
295 raise RelaxLenError('increment', n)
296 elif isinstance(inc, int):
297 inc = [inc]*n
298
299
300 if not lower:
301
302 lower = []
303 upper = []
304
305
306 if cdp.model_type == 'R2eff':
307
308 for exp_type, frq, offset, point in loop_exp_frq_offset_point():
309
310 for param_name, param_index, si, r20_key in loop_parameters(spins=spins):
311
312 if param_name == 'r2eff':
313 lower.append(1.0)
314 upper.append(40.0)
315
316
317 elif param_name == 'i0':
318 lower.append(0.0001)
319 upper.append(max(spins[si].peak_intensity.values()))
320
321
322 else:
323
324 for param_name, param_index, si, r20_key in loop_parameters(spins=spins):
325
326 if si == None:
327 si = 0
328
329
330 if param_name in ['r2', 'r2a', 'r2b']:
331 lower.append(5.0)
332 upper.append(20.0)
333
334
335 elif param_name in ['phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2']:
336 lower.append(0.0)
337 upper.append(10.0)
338
339
340 elif param_name in ['dw', 'dw_AB', 'dw_AC', 'dw_BC']:
341 if spins[si].model in MODEL_LIST_MMQ:
342 lower.append(-10.0)
343 else:
344 lower.append(0.0)
345 upper.append(10.0)
346
347
348 elif param_name in ['dwH', 'dwH_AB', 'dwH_AC', 'dwH_BC']:
349 if spins[si].model in MODEL_LIST_MMQ:
350 lower.append(-3.0)
351 else:
352 lower.append(0.0)
353 upper.append(3.0)
354
355
356 elif param_name == 'pA':
357 if spins[si].model == MODEL_M61B:
358 lower.append(0.85)
359 else:
360 lower.append(0.5)
361 upper.append(1.0)
362
363
364 elif param_name == 'pB':
365 lower.append(0.0)
366 upper.append(0.5)
367
368
369 elif param_name in ['kex', 'kex_AB', 'kex_AC', 'kex_BC', 'k_AB', 'kB', 'kC']:
370 lower.append(1.0)
371 upper.append(10000.0)
372
373
374 elif param_name in ['tex']:
375 lower.append(1/10000.0)
376 upper.append(1.0)
377
378
379 for param_name, param_index, si, r20_key in loop_parameters(spins=spins):
380
381 if si == None:
382 si = 0
383
384
385 if hasattr(spins[si], param_name):
386 val = getattr(spins[si], param_name)
387
388
389 if is_float(val) and val != 0.0:
390
391 print("The spin '%s' parameter '%s' is pre-set to %s, skipping it in the grid search." % (spin_ids[si], param_name, val))
392
393
394 inc[param_index] = 1
395 lower[param_index] = val
396 upper[param_index] = val
397
398
399 if isinstance(val, dict) and len(val) > 0:
400
401 if r20_key != None:
402 try:
403 val_dic = val[r20_key]
404 except KeyError:
405 print("The key:%s does not exist"%r20_key)
406 continue
407
408 if is_float(val_dic):
409
410 print("The spin '%s' parameter %s '%s[%i]' is pre-set to %s, skipping it in the grid search." % (spin_ids[si], r20_key, param_name, param_index, val_dic))
411
412
413 inc[param_index] = 1
414 lower[param_index] = val_dic
415 upper[param_index] = val_dic
416
417
418 grid_size = 1
419 for i in range(n):
420 grid_size *= inc[i]
421
422
423 lower_new = []
424 upper_new = []
425 for i in range(n):
426 lower_new.append(lower[i] / scaling_matrix[i, i])
427 upper_new.append(upper[i] / scaling_matrix[i, i])
428
429
430 return grid_size, inc, lower_new, upper_new
431
432
433 -def minimise_r2eff(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):
434 """Optimise the R2eff model by fitting the 2-parameter exponential curves.
435
436 This mimics the R1 and R2 relax_fit analysis.
437
438
439 @keyword min_algor: The minimisation algorithm to use.
440 @type min_algor: str
441 @keyword min_options: An array of options to be used by the minimisation algorithm.
442 @type min_options: array of str
443 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
444 @type func_tol: None or float
445 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
446 @type grad_tol: None or float
447 @keyword max_iterations: The maximum number of iterations for the algorithm.
448 @type max_iterations: int
449 @keyword constraints: If True, constraints are used during optimisation.
450 @type constraints: bool
451 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned.
452 @type scaling: bool
453 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
454 @type verbosity: int
455 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
456 @type sim_index: None or int
457 @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.
458 @type lower: array of numbers
459 @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.
460 @type upper: array of numbers
461 @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.
462 @type inc: array of int
463 """
464
465
466 if not C_module_exp_fn:
467 raise RelaxError("Relaxation curve fitting is not available. Try compiling the C modules on your platform.")
468
469
470 for spin, spin_id in spin_loop(return_id=True, skip_desel=True):
471
472 if not hasattr(spin, 'peak_intensity'):
473 continue
474
475
476 for exp_type, frq, offset, point in loop_exp_frq_offset_point():
477
478 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point)
479
480
481 param_vector = assemble_param_vector(spins=[spin], key=param_key, sim_index=sim_index)
482
483
484 scaling_matrix = assemble_scaling_matrix(spins=[spin], key=param_key, scaling=scaling)
485 if len(scaling_matrix):
486 param_vector = dot(inv(scaling_matrix), param_vector)
487
488
489 lower_new, upper_new = None, None
490 if match('^[Gg]rid', min_algor):
491 grid_size, inc_new, lower_new, upper_new = grid_search_setup(spins=[spin], spin_ids=[spin_id], param_vector=param_vector, lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix)
492
493
494 A, b = None, None
495 if constraints:
496 A, b = linear_constraints(spins=[spin], scaling_matrix=scaling_matrix)
497
498
499 if verbosity >= 1:
500
501 top = 2
502 if verbosity >= 2:
503 top += 2
504 text = "Fitting to spin %s, frequency %s and dispersion point %s" % (spin_id, frq, point)
505 subsection(file=sys.stdout, text=text, prespace=top)
506
507
508 if match('^[Gg]rid', min_algor):
509 print("Unconstrained grid search size: %s (constraints may decrease this size).\n" % grid_size)
510
511
512 values = []
513 errors = []
514 times = []
515 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point):
516 values.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, sim_index=sim_index))
517 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True))
518 times.append(time)
519
520
521 scaling_list = []
522 for i in range(len(scaling_matrix)):
523 scaling_list.append(scaling_matrix[i, i])
524
525
526 setup(num_params=len(param_vector), num_times=len(times), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list)
527
528
529 if search('^[Gg]rid', min_algor):
530 results = grid(func=func, args=(), num_incs=inc_new, lower=lower_new, upper=upper_new, A=A, b=b, verbosity=verbosity)
531
532
533 param_vector, chi2, iter_count, warning = results
534 f_count = iter_count
535 g_count = 0.0
536 h_count = 0.0
537
538
539 else:
540 results = generic_minimise(func=func, dfunc=dfunc, d2func=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)
541
542
543 if results == None:
544 return
545 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results
546
547
548 if scaling:
549 param_vector = dot(scaling_matrix, param_vector)
550
551
552 disassemble_param_vector(param_vector=param_vector, spins=[spin], key=param_key, sim_index=sim_index)
553
554
555 if sim_index != None:
556
557 spin.chi2_sim[sim_index] = chi2
558
559
560 spin.iter_sim[sim_index] = iter_count
561
562
563 spin.f_count_sim[sim_index] = f_count
564
565
566 spin.g_count_sim[sim_index] = g_count
567
568
569 spin.h_count_sim[sim_index] = h_count
570
571
572 spin.warning_sim[sim_index] = warning
573
574
575 else:
576
577 spin.chi2 = chi2
578
579
580 spin.iter = iter_count
581
582
583 spin.f_count = f_count
584
585
586 spin.g_count = g_count
587
588
589 spin.h_count = h_count
590
591
592 spin.warning = warning
593
594
595
596
598 """The relaxation dispersion memo class."""
599
600 - def __init__(self, spins=None, spin_ids=None, sim_index=None, scaling_matrix=None, verbosity=None):
601 """Initialise the relaxation dispersion memo class.
602
603 This is used for handling the optimisation results returned from a slave processor. It runs on the master processor and is used to store data which is passed to the slave processor and then passed back to the master via the results command.
604
605
606 @keyword spins: The list of spin data container for the cluster. If this argument is supplied, then the spin_id argument will be ignored.
607 @type spins: list of SpinContainer instances
608 @keyword spin_ids: The spin ID strings for the cluster.
609 @type spin_ids: list of str
610 @keyword sim_index: The optional MC simulation index.
611 @type sim_index: int
612 @keyword scaling_matrix: The diagonal, square scaling matrix.
613 @type scaling_matrix: numpy diagonal matrix
614 @keyword verbosity: The verbosity level. This is used by the result command returned to the master for printouts.
615 @type verbosity: int
616 """
617
618
619 super(Disp_memo, self).__init__()
620
621
622 self.spins = spins
623 self.spin_ids = spin_ids
624 self.sim_index = sim_index
625 self.scaling_matrix = scaling_matrix
626 self.verbosity = verbosity
627
628
629
631 """Command class for relaxation dispersion optimisation on the slave processor."""
632
633 - def __init__(self, spins=None, spin_ids=None, sim_index=None, scaling_matrix=None, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, verbosity=0, lower=None, upper=None, inc=None, fields=None, param_names=None):
634 """Initialise the base class, storing all the master data to be sent to the slave processor.
635
636 This method is run on the master processor whereas the run() method is run on the slave processor.
637
638
639 @keyword spins: The list of spin data container for the cluster. If this argument is supplied, then the spin_id argument will be ignored.
640 @type spins: list of SpinContainer instances
641 @keyword spin_ids: The list of spin ID strings corresponding to the spins argument.
642 @type spin_ids: list of str
643 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
644 @type sim_index: None or int
645 @keyword scaling_matrix: The diagonal, square scaling matrix.
646 @type scaling_matrix: numpy diagonal matrix
647 @keyword min_algor: The minimisation algorithm to use.
648 @type min_algor: str
649 @keyword min_options: An array of options to be used by the minimisation algorithm.
650 @type min_options: array of str
651 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
652 @type func_tol: None or float
653 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
654 @type grad_tol: None or float
655 @keyword max_iterations: The maximum number of iterations for the algorithm.
656 @type max_iterations: int
657 @keyword constraints: If True, constraints are used during optimisation.
658 @type constraints: bool
659 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
660 @type verbosity: int
661 @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.
662 @type lower: array of numbers
663 @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.
664 @type upper: array of numbers
665 @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.
666 @type inc: array of int
667 @keyword fields: The list of unique of spectrometer field strengths.
668 @type fields: int
669 @keyword param_names: The list of parameter names to use in printouts.
670 @type param_names: str
671 """
672
673
674 super(Disp_minimise_command, self).__init__()
675
676
677 self.spins = spins
678 self.spin_ids = spin_ids
679 self.sim_index = sim_index
680 self.scaling_matrix = scaling_matrix
681 self.verbosity = verbosity
682 self.min_algor = min_algor
683 self.min_options = min_options
684 self.func_tol = func_tol
685 self.grad_tol = grad_tol
686 self.max_iterations = max_iterations
687 self.fields = fields
688 self.param_names = param_names
689
690
691 self.param_vector = assemble_param_vector(spins=self.spins)
692 if len(scaling_matrix):
693 self.param_vector = dot(inv(scaling_matrix), self.param_vector)
694
695
696 self.lower_new, self.upper_new = None, None
697 if search('^[Gg]rid', min_algor):
698 self.grid_size, self.inc_new, self.lower_new, self.upper_new = grid_search_setup(spins=spins, spin_ids=spin_ids, param_vector=self.param_vector, lower=lower, upper=upper, inc=inc, scaling_matrix=self.scaling_matrix)
699
700
701 self.A, self.b = None, None
702 if constraints:
703 self.A, self.b = linear_constraints(spins=spins, scaling_matrix=scaling_matrix)
704
705
706 if spins[0].model in [MODEL_LM63, MODEL_CR72, MODEL_CR72_FULL, MODEL_M61, MODEL_TP02, MODEL_TAP03, MODEL_MP05] and not hasattr(cdp, 'spectrometer_frq'):
707 raise RelaxError("The spectrometer frequency information has not been specified.")
708
709
710 self.values, self.errors, self.missing, self.frqs, self.frqs_H, self.exp_types, self.relax_times = return_r2eff_arrays(spins=spins, spin_ids=spin_ids, fields=fields, field_count=len(fields), sim_index=sim_index)
711
712
713 self.chemical_shifts, self.offsets, self.tilt_angles, self.Delta_omega, self.w_eff = return_offset_data(spins=spins, spin_ids=spin_ids, field_count=len(fields))
714 self.r1 = return_r1_data(spins=spins, spin_ids=spin_ids, field_count=len(fields), sim_index=sim_index)
715
716
717 self.param_num = param_num(spins=spins)
718
719
720 self.dispersion_points = cdp.dispersion_points
721 self.cpmg_frqs = return_cpmg_frqs(ref_flag=False)
722 self.spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False)
723
724
725 - def run(self, processor, completed):
726 """Set up and perform the optimisation."""
727
728
729 if self.verbosity >= 1:
730
731 top = 2
732 if self.verbosity >= 2:
733 top += 2
734 subsection(file=sys.stdout, text="Fitting to the spin block %s"%self.spin_ids, prespace=top)
735
736
737 if search('^[Gg]rid', self.min_algor):
738 print("Unconstrained grid search size: %s (constraints may decrease this size).\n" % self.grid_size)
739
740
741 model = Dispersion(model=self.spins[0].model, num_params=self.param_num, num_spins=count_spins(self.spins), num_frq=len(self.fields), exp_types=self.exp_types, values=self.values, errors=self.errors, missing=self.missing, frqs=self.frqs, frqs_H=self.frqs_H, cpmg_frqs=self.cpmg_frqs, spin_lock_nu1=self.spin_lock_nu1, chemical_shifts=self.chemical_shifts, offset=self.offsets, tilt_angles=self.tilt_angles, r1=self.r1, relax_times=self.relax_times, scaling_matrix=self.scaling_matrix)
742
743
744 if search('^[Gg]rid', self.min_algor):
745 results = grid(func=model.func, args=(), num_incs=self.inc_new, lower=self.lower_new, upper=self.upper_new, A=self.A, b=self.b, verbosity=self.verbosity)
746
747
748 param_vector, chi2, iter_count, warning = results
749 f_count = iter_count
750 g_count = 0.0
751 h_count = 0.0
752
753
754 else:
755 results = generic_minimise(func=model.func, args=(), x0=self.param_vector, min_algor=self.min_algor, min_options=self.min_options, func_tol=self.func_tol, grad_tol=self.grad_tol, maxiter=self.max_iterations, A=self.A, b=self.b, full_output=True, print_flag=self.verbosity)
756
757
758 if results == None:
759 return
760 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results
761
762
763 if self.verbosity:
764 print("\nOptimised parameter values:")
765 for i in range(len(param_vector)):
766 print("%-20s %25.15f" % (self.param_names[i], param_vector[i]*self.scaling_matrix[i, i]))
767
768
769 if self.sim_index != None:
770 print("Simulation %s, cluster %s" % (self.sim_index+1, self.spin_ids))
771
772
773 processor.return_object(Disp_result_command(processor=processor, memo_id=self.memo_id, param_vector=param_vector, chi2=chi2, iter_count=iter_count, f_count=f_count, g_count=g_count, h_count=h_count, warning=warning, missing=self.missing, back_calc=model.back_calc, completed=False))
774
775
776
778 """Class for processing the dispersion optimisation results.
779
780 This object will be sent from the slave back to the master to have its run() method executed.
781 """
782
783 - def __init__(self, processor=None, memo_id=None, param_vector=None, chi2=None, iter_count=None, f_count=None, g_count=None, h_count=None, warning=None, missing=None, back_calc=None, completed=True):
784 """Set up this class object on the slave, placing the minimisation results here.
785
786 @keyword processor: The processor object.
787 @type processor: multi.processor.Processor instance
788 @keyword memo_id: The memo identification string.
789 @type memo_id: str
790 @keyword param_vector: The optimised parameter vector.
791 @type param_vector: numpy rank-1 array
792 @keyword chi2: The final target function value.
793 @type chi2: float
794 @keyword iter_count: The number of optimisation iterations.
795 @type iter_count: int
796 @keyword f_count: The total function call count.
797 @type f_count: int
798 @keyword g_count: The total gradient call count.
799 @type g_count: int
800 @keyword h_count: The total Hessian call count.
801 @type h_count: int
802 @keyword warning: Any optimisation warnings.
803 @type warning: str or None
804 @keyword missing: The data structure indicating which R2eff/R1rho' base data is missing.
805 @type missing: numpy rank-3 array
806 @keyword back_calc: The back-calculated R2eff/R1rho' data structure from the target function class. This is will be transfered to the master to be stored in the r2eff_bc data structure.
807 @type back_calc: numpy rank-3 array
808 @keyword completed: A flag which if True signals that the optimisation successfully completed.
809 @type completed: bool
810 """
811
812
813 super(Disp_result_command, self).__init__(processor=processor, completed=completed)
814
815
816 self.memo_id = memo_id
817 self.param_vector = param_vector
818 self.chi2 = chi2
819 self.iter_count = iter_count
820 self.f_count = f_count
821 self.g_count = g_count
822 self.h_count = h_count
823 self.warning = warning
824 self.missing = missing
825 self.back_calc = back_calc
826 self.completed = completed
827
828
829 - def run(self, processor=None, memo=None):
830 """Disassemble the model-free optimisation results (on the master).
831
832 @param processor: Unused!
833 @type processor: None
834 @param memo: The dispersion memo. This holds a lot of the data and objects needed for processing the results from the slave.
835 @type memo: memo
836 """
837
838
839 if memo.scaling_matrix != None:
840 param_vector = dot(memo.scaling_matrix, self.param_vector)
841
842
843 disassemble_param_vector(param_vector=param_vector, spins=memo.spins, sim_index=memo.sim_index)
844 param_conversion(spins=memo.spins, sim_index=memo.sim_index)
845
846
847 if memo.sim_index != None:
848 for spin in memo.spins:
849
850 if not spin.select:
851 continue
852
853
854 spin.chi2_sim[memo.sim_index] = self.chi2
855
856
857 spin.iter_sim[memo.sim_index] = self.iter_count
858
859
860 spin.f_count_sim[memo.sim_index] = self.f_count
861
862
863 spin.g_count_sim[memo.sim_index] = self.g_count
864
865
866 spin.h_count_sim[memo.sim_index] = self.h_count
867
868
869 spin.warning_sim[memo.sim_index] = self.warning
870
871
872 else:
873 for spin in memo.spins:
874
875 if not spin.select:
876 continue
877
878
879 spin.chi2 = self.chi2
880
881
882 spin.iter = self.iter_count
883
884
885 spin.f_count = self.f_count
886
887
888 spin.g_count = self.g_count
889
890
891 spin.h_count = self.h_count
892
893
894 spin.warning = self.warning
895
896
897 if memo.sim_index == None:
898
899 proton_mmq_flag = has_proton_mmq_cpmg()
900
901
902 si = 0
903 for spin_index in range(len(memo.spins)):
904
905 if not memo.spins[spin_index].select:
906 continue
907
908
909 pack_back_calc_r2eff(spin=memo.spins[spin_index], spin_id=memo.spin_ids[spin_index], si=si, back_calc=self.back_calc, proton_mmq_flag=proton_mmq_flag)
910
911
912 si += 1
913