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 relaxation dispersion API object."""
26
27
28 import bmrblib
29 from copy import deepcopy
30 from numpy import int32, sqrt, zeros
31 from re import match, search
32 import string
33 import sys
34 from types import MethodType
35
36
37 from lib.arg_check import is_list, is_str_list
38 from lib.dispersion.variables import EXP_TYPE_CPMG_PROTON_MQ, EXP_TYPE_CPMG_PROTON_SQ, MODEL_LIST_MMQ, MODEL_R2EFF, PARAMS_R20
39 from lib.errors import RelaxError, RelaxImplementError
40 from lib.text.sectioning import subsection
41 from multi import Processor_box
42 from pipe_control import pipes, sequence
43 from pipe_control.exp_info import bmrb_write_citations, bmrb_write_methods, bmrb_write_software
44 from pipe_control.mol_res_spin import bmrb_write_entity, check_mol_res_spin_data, get_molecule_names, return_spin, spin_loop
45 from pipe_control.pipes import check_pipe
46 from pipe_control.spectrometer import check_spectrometer_setup
47 from pipe_control.sequence import return_attached_protons
48 from specific_analyses.api_base import API_base
49 from specific_analyses.api_common import API_common
50 from specific_analyses.relax_disp.checks import check_model_type
51 from specific_analyses.relax_disp.data import average_intensity, calc_rotating_frame_params, find_intensity_keys, generate_r20_key, has_exponential_exp_type, has_proton_mmq_cpmg, loop_cluster, loop_exp_frq, loop_exp_frq_offset_point, loop_time, pack_back_calc_r2eff, return_param_key_from_data, spin_ids_to_containers
52 from specific_analyses.relax_disp.optimisation import Disp_memo, Disp_minimise_command, back_calc_peak_intensities, back_calc_r2eff, calculate_r2eff, minimise_r2eff
53 from specific_analyses.relax_disp.parameter_object import Relax_disp_params
54 from specific_analyses.relax_disp.parameters import get_param_names, get_value, loop_parameters, param_index_to_param_info, param_num, r1_setup
55
56
58 """Class containing functions for relaxation dispersion curve fitting."""
59
60
61 instance = None
62
74
75
77 """Custom generator method for looping over the base data.
78
79 For the R2eff model, the base data is the peak intensity data defining a single exponential curve. For all other models, the base data is the R2eff/R1rho values for individual spins.
80
81
82 @return: For the R2eff model, a tuple of the spin container and the exponential curve identifying key (the CPMG frequency or R1rho spin-lock field strength). For all other models, the spin container and ID from the spin loop.
83 @rtype: (tuple of SpinContainer instance and float) or (SpinContainer instance and str)
84 """
85
86
87 if cdp.model_type == MODEL_R2EFF:
88
89 for spin, spin_id in spin_loop(return_id=True):
90
91 if not spin.select:
92 continue
93
94
95 if not hasattr(spin, 'peak_intensity'):
96 continue
97
98
99 for exp_type, frq, offset, point in loop_exp_frq_offset_point():
100 yield spin, spin_id, exp_type, frq, offset, point
101
102
103 else:
104
105 proton_mmq_flag = has_proton_mmq_cpmg()
106
107
108 for spin, spin_id in spin_loop(return_id=True):
109
110 if not spin.select:
111 continue
112
113
114 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H':
115 continue
116
117
118 proton = None
119 if proton_mmq_flag:
120 proton = return_attached_protons(spin_id)[0]
121
122
123 if not hasattr(spin, 'r2eff') and not hasattr(proton, 'r2eff'):
124 continue
125
126
127 yield spin, spin_id
128
129
131 """Write the model-free results to a BMRB NMR-STAR v3.1 formatted file.
132
133 @param file_path: The full file path.
134 @type file_path: str
135 @keyword version: The BMRB NMR-STAR dictionary format to output to.
136 @type version: str
137 """
138
139
140
141
142
143 raise RelaxImplementError('bmrb_write')
144
145
146 check_spectrometer_setup(escalate=2)
147
148
149 cdp = pipes.get_pipe()
150
151
152 star = bmrblib.create_nmr_star('relax_relaxation_dispersion_results', file_path, version)
153
154
155 mol_name_list = []
156 res_num_list = []
157 res_name_list = []
158 atom_name_list = []
159
160 isotope_list = []
161 element_list = []
162
163 chi2_list = []
164 model_list = []
165
166
167 for spin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True):
168
169 if spin.name == 'H' or (hasattr(spin, 'element') and spin.element == 'H'):
170 warn(RelaxWarning("Skipping the proton spin '%s'." % spin_id))
171 continue
172
173
174 if res_num == None:
175 raise RelaxError("For the BMRB, the residue of spin '%s' must be numbered." % spin_id)
176 if res_name == None:
177 raise RelaxError("For the BMRB, the residue of spin '%s' must be named." % spin_id)
178 if spin.name == None:
179 raise RelaxError("For the BMRB, the spin '%s' must be named." % spin_id)
180 if not hasattr(spin, 'isotope') or spin.isotope == None:
181 raise RelaxError("For the BMRB, the spin isotope type of '%s' must be specified." % spin_id)
182 if not hasattr(spin, 'element') or spin.element == None:
183 raise RelaxError("For the BMRB, the spin element type of '%s' must be specified. Please use the spin user function for setting the element type." % spin_id)
184
185
186 mol_name_list.append(mol_name)
187 res_num_list.append(res_num)
188 res_name_list.append(res_name)
189 atom_name_list.append(spin.name)
190
191
192 if hasattr(spin, 'isotope'):
193 isotope_list.append(int(spin.isotope.strip(string.ascii_letters)))
194 else:
195 isotope_list.append(None)
196
197
198 if hasattr(spin, 'element'):
199 element_list.append(spin.element)
200 else:
201 element_list.append(None)
202
203
204 if hasattr(spin, 'chi2'):
205 chi2_list.append(spin.chi2)
206 else:
207 chi2_list.append(None)
208
209
210 model_list.append(spin.model)
211
212
213 entity_ids = zeros(len(mol_name_list), int32)
214 mol_names = get_molecule_names()
215 for i in range(len(mol_name_list)):
216 for j in range(len(mol_names)):
217 if mol_name_list[i] == mol_names[j]:
218 entity_ids[i] = j+1
219
220
221
222
223
224
225 bmrb_write_citations(star)
226
227
228
229
230
231
232 bmrb_write_entity(star)
233
234
235
236
237
238
239 bmrb_write_methods(star)
240
241
242 software_ids, software_labels = bmrb_write_software(star)
243
244
245 - def calculate(self, spin_id=None, scaling_matrix=None, verbosity=1, sim_index=None):
246 """Calculate the model chi-squared value or the R2eff values for fixed time period data.
247
248 @keyword spin_id: The spin identification string.
249 @type spin_id: None or str
250 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices.
251 @type scaling_matrix: list of numpy rank-2, float64 array or list of None
252 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
253 @type verbosity: int
254 @keyword sim_index: The MC simulation index (unused).
255 @type sim_index: None
256 """
257
258
259 check_pipe()
260 check_mol_res_spin_data()
261 check_model_type()
262
263
264 if cdp.model_type == MODEL_R2EFF:
265 calculate_r2eff()
266
267
268 else:
269
270 proton_mmq_flag = has_proton_mmq_cpmg()
271
272
273 model_index = -1
274 for spin_ids in self.model_loop():
275
276 model_index += 1
277
278
279 spins = spin_ids_to_containers(spin_ids)
280
281
282 skip = True
283 for spin in spins:
284 if spin.select:
285 skip = False
286 if skip:
287 continue
288
289
290 back_calc = back_calc_r2eff(spins=spins, spin_ids=spin_ids, store_chi2=True)
291
292
293 for i, spin_id in enumerate(spin_ids):
294 spin = spins[i]
295 pack_back_calc_r2eff(spin=spin, spin_id=spin_id, si=i, back_calc=back_calc, proton_mmq_flag=proton_mmq_flag)
296
297
299 """Return the 'Log barrier' optimisation constraint algorithm.
300
301 @return: The 'Log barrier' constraint algorithm.
302 @rtype: str
303 """
304
305
306 return 'Log barrier'
307
308
310 """Create the Monte Carlo peak intensity data.
311
312 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method.
313 @type data_id: SpinContainer instance and float
314 @return: The Monte Carlo simulation data.
315 @rtype: list of floats
316 """
317
318
319 if cdp.model_type == MODEL_R2EFF:
320
321 spin, spin_id, exp_type, frq, offset, point = data_id
322
323
324 values = back_calc_peak_intensities(spin=spin, spin_id=spin_id, exp_type=exp_type, frq=frq, offset=offset, point=point)
325
326
327 else:
328
329 proton_mmq_flag = has_proton_mmq_cpmg()
330
331
332 spin, spin_id = data_id
333
334
335 back_calc = back_calc_r2eff(spins=[spin], spin_ids=[spin_id])
336
337
338 if proton_mmq_flag:
339 proton = return_attached_protons(spin_id)[0]
340 proton_back_calc = back_calc_r2eff(spins=[proton], spin_ids=[spin_id])
341
342
343 values = {}
344 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True):
345
346 current_bc = back_calc
347 current_spin = spin
348 if exp_type in [EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_PROTON_MQ]:
349 current_spin = proton
350 current_bc = proton_back_calc
351
352
353 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point)
354
355
356 if not hasattr(current_spin, 'r2eff') or param_key not in current_spin.r2eff:
357 continue
358
359
360 values[param_key] = back_calc[ei][0][mi][oi][di]
361
362
363 return values
364
365
366 - def deselect(self, sim_index=None, model_info=None):
367 """Deselect models or simulations.
368
369 @keyword sim_index: The optional Monte Carlo simulation index. If None, then models will be deselected, otherwise the given simulation will.
370 @type sim_index: None or int
371 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
372 @type model_info: list of SpinContainer instances, list of str
373 """
374
375
376 for spin_id in model_info:
377
378 spin = return_spin(spin_id)
379
380
381 if sim_index == None:
382 spin.select = False
383
384
385 else:
386 spin.select_sim[sim_index] = False
387
388
389 - def duplicate_data(self, pipe_from=None, pipe_to=None, model_info=None, global_stats=False, verbose=True):
390 """Duplicate the data specific to a single model.
391
392 @keyword pipe_from: The data pipe to copy the data from.
393 @type pipe_from: str
394 @keyword pipe_to: The data pipe to copy the data to.
395 @type pipe_to: str
396 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
397 @type model_info: list of SpinContainer instances, list of str
398 @keyword global_stats: The global statistics flag.
399 @type global_stats: bool
400 @keyword verbose: A flag which if True will cause info to be printed out.
401 @type verbose: bool
402 """
403
404
405 if not pipes.has_pipe(pipe_to):
406 pipes.create(pipe_to, pipe_type='relax_disp', switch=False)
407
408
409 dp_from = pipes.get_pipe(pipe_from)
410 dp_to = pipes.get_pipe(pipe_to)
411
412
413 for data_name in dir(dp_from):
414
415 if data_name in ['mol', 'interatomic', 'structure', 'exp_info', 'result_files']:
416 continue
417
418
419 if data_name in ['model']:
420 continue
421
422
423 if search('^_', data_name) or data_name in dp_from.__class__.__dict__:
424 continue
425
426
427 data_from = getattr(dp_from, data_name)
428
429
430 if hasattr(dp_to, data_name):
431
432 data_to = getattr(dp_to, data_name)
433
434
435 if data_from != data_to:
436 raise RelaxError("The object " + repr(data_name) + " is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
437
438
439 continue
440
441
442 setattr(dp_to, data_name, deepcopy(data_from))
443
444
445 if dp_from.mol.is_empty():
446 return
447
448
449 if dp_to.mol.is_empty():
450 sequence.copy(pipe_from=pipe_from, pipe_to=pipe_to, preserve_select=True, verbose=verbose)
451
452
453 for id in model_info:
454
455 spin = return_spin(id, pipe=pipe_from)
456
457
458 for name in dir(spin):
459
460 if search('^__', name):
461 continue
462
463
464 obj = getattr(spin, name)
465
466
467 if isinstance(obj, MethodType):
468 continue
469
470
471 new_obj = deepcopy(getattr(spin, name))
472 setattr(dp_to.mol[spin._mol_index].res[spin._res_index].spin[spin._spin_index], name, new_obj)
473
474
475 - def eliminate(self, name, value, args, sim=None, model_info=None):
476 """Relaxation dispersion model elimination, parameter by parameter.
477
478 @param name: The parameter name.
479 @type name: str
480 @param value: The parameter value.
481 @type value: float
482 @param args: The c1 and c2 elimination constant overrides.
483 @type args: None or tuple of float
484 @keyword sim: The Monte Carlo simulation index.
485 @type sim: int
486 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
487 @type model_info: list of SpinContainer instances, list of str
488 @return: True if the model is to be eliminated, False otherwise.
489 @rtype: bool
490 """
491
492
493 if name in ['r2eff', 'i0']:
494 return False
495
496
497 c1 = 0.501
498 c2 = 0.999
499 c3 = 1.0
500
501
502 if args != None:
503 c1, c2, c3 = args
504
505
506 elim_text = "Data pipe '%s': The %s parameter of %.5f is %s, eliminating " % (pipes.cdp_name(), name, value, "%s")
507 if sim == None:
508 elim_text += "the spin cluster %s." % model_info
509 else:
510 elim_text += "simulation %i of the spin cluster %s." % (sim, model_info)
511
512
513 if name == 'pA':
514 if value < c1:
515 print(elim_text % ("less than %.5f" % c1))
516 return True
517 if value > c2:
518 print(elim_text % ("greater than %.5f" % c2))
519 return True
520
521
522 if name == 'tex':
523 if value > c3:
524 print(elim_text % ("greater than %.5f" % c3))
525 return True
526
527
528 return False
529
530
532 """Return a vector of parameter names.
533
534 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
535 @type model_info: list of SpinContainer instances, list of str
536 @return: The vector of parameter names.
537 @rtype: list of str
538 """
539
540
541 spins = []
542 for spin_id in model_info:
543
544 spin = return_spin(spin_id)
545
546
547 if not spin.select:
548 continue
549
550
551 spins.append(spin)
552
553
554 if not len(spins):
555 return None
556
557
558 return get_param_names(spins)
559
560
562 """Return a vector of parameter values.
563
564 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
565 @type model_info: list of SpinContainer instances, list of str
566 @keyword sim_index: The Monte Carlo simulation index.
567 @type sim_index: int
568 @return: The vector of parameter values.
569 @rtype: list of str
570 """
571
572
573 spins = []
574 for spin_id in model_info:
575
576 spin = return_spin(spin_id)
577
578
579 if not spin.select:
580 continue
581
582
583 spins.append(spin)
584
585
586 if not len(spins):
587 return None
588
589
590 r1_setup()
591
592
593 values = []
594 for param_name, param_index, si, r20_key in loop_parameters(spins=spins):
595 values.append(get_value(spins=spins, sim_index=sim_index, param_name=param_name, spin_index=si, r20_key=r20_key))
596
597
598 return values
599
600
601 - def grid_search(self, lower=None, upper=None, inc=None, scaling_matrix=None, constraints=True, verbosity=1, sim_index=None):
602 """The relaxation dispersion curve fitting grid search function.
603
604 @keyword lower: The per-model lower bounds of the grid search which must be equal to the number of parameters in the model.
605 @type lower: list of list of numbers
606 @keyword upper: The per-model upper bounds of the grid search which must be equal to the number of parameters in the model.
607 @type upper: list of list of numbers
608 @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.
609 @type inc: list of list of int
610 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices.
611 @type scaling_matrix: list of numpy rank-2, float64 array or list of None
612 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used.
613 @type constraints: bool
614 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity.
615 @type verbosity: int
616 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised.
617 @type sim_index: int
618 """
619
620
621 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
622
623
625 """Create bounds for the OpenDX mapping function.
626
627 @param param: The name of the parameter to return the lower and upper bounds of.
628 @type param: str
629 @param spin_id: The spin identification string (unused).
630 @type spin_id: None
631 @return: The upper and lower bounds of the parameter.
632 @rtype: list of float
633 """
634
635
636 if not self._PARAMS.contains(param):
637 raise RelaxError("The parameter '%s' is not valid for this data pipe type." % param)
638
639
640 spin = return_spin(spin_id)
641
642
643 param_keys = []
644 for exp_type, frq, offset, point in loop_exp_frq_offset_point():
645
646 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point)
647
648
649 param_keys.append(param_key)
650
651
652 param_vector = []
653
654
655 param_names = []
656 for param_name, param_index, si, r20_key in loop_parameters(spins=[spin]):
657
658 param_vector.append([0.0])
659
660
661 param_names.append(param_name)
662
663
664 for i in range(len(param_names)):
665
666
667 if param_names[i] == param:
668 return [self._PARAMS.grid_lower(param, incs=0, model_info=[spin_id]), self._PARAMS.grid_upper(param, incs=0, model_info=[spin_id])]
669
670
671 - 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):
672 """Relaxation dispersion curve fitting function.
673
674 @keyword min_algor: The minimisation algorithm to use.
675 @type min_algor: str
676 @keyword min_options: An array of options to be used by the minimisation algorithm.
677 @type min_options: array of str
678 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
679 @type func_tol: None or float
680 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
681 @type grad_tol: None or float
682 @keyword max_iterations: The maximum number of iterations for the algorithm.
683 @type max_iterations: int
684 @keyword constraints: If True, constraints are used during optimisation.
685 @type constraints: bool
686 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices.
687 @type scaling_matrix: list of numpy rank-2, float64 array or list of None
688 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
689 @type verbosity: int
690 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
691 @type sim_index: None or int
692 @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.
693 @type lower: list of list of numbers
694 @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.
695 @type upper: list of list of numbers
696 @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.
697 @type inc: list of list of int
698 """
699
700
701 check_mol_res_spin_data()
702 check_model_type()
703
704
705 algor = min_algor
706 if min_algor == 'Log barrier':
707 algor = min_options[0]
708
709 allow = False
710
711 if hasattr(cdp, 'model_type'):
712
713 model_type = cdp.model_type
714
715 if model_type == MODEL_R2EFF:
716 if match('^[Gg]rid$', algor):
717 allow = True
718
719 elif match('^[Ss]implex$', algor):
720 allow = True
721
722
723 elif match('^[Bb][Ff][Gg][Ss]$', algor):
724 allow = True
725
726
727 elif match('^[Nn]ewton$', algor):
728 allow = True
729
730
731 elif match('^[Nn]ewton$', algor):
732 allow = True
733
734
735 elif match('^[Mm][Oo][Mm]$', algor) or match('[Mm]ethod of [Mm]ultipliers$', algor):
736 allow = True
737
738
739 elif match('^[Ll]og [Bb]arrier$', algor):
740 allow = True
741
742
743 else:
744 if match('^[Gg]rid$', algor):
745 allow = True
746
747 elif match('^[Ss]implex$', algor):
748 allow = True
749
750
751 else:
752 model_type = 'None'
753
754 allow = False
755
756 if not allow:
757 raise RelaxError("Minimisation algorithm '%s' is not allowed, since function gradients for model '%s' is not implemented. Only the 'simplex' minimisation algorithm is supported for the relaxation dispersion analysis of this model."%(algor, model_type))
758
759
760 if not hasattr(cdp, 'cpmg_frqs_list'):
761 cdp.cpmg_frqs_list = []
762 if not hasattr(cdp, 'spin_lock_nu1_list'):
763 cdp.spin_lock_nu1_list = []
764
765
766 processor_box = Processor_box()
767 processor = processor_box.processor
768
769
770 num_time_pts = 1
771 if hasattr(cdp, 'num_time_pts'):
772 num_time_pts = cdp.num_time_pts
773
774
775 fields = [None]
776 field_count = 1
777 if hasattr(cdp, 'spectrometer_frq'):
778 fields = cdp.spectrometer_frq_list
779 field_count = cdp.spectrometer_frq_count
780
781
782 model_index = -1
783 for spin_ids in self.model_loop():
784
785 model_index += 1
786
787
788 spins = spin_ids_to_containers(spin_ids)
789
790
791 skip = True
792 for spin in spins:
793 if spin.select:
794 skip = False
795 if skip:
796 continue
797
798
799 lower_i, upper_i, inc_i = None, None, None
800 if min_algor == 'grid':
801 lower_i = lower[model_index]
802 upper_i = upper[model_index]
803 inc_i = inc[model_index]
804
805
806 if cdp.model_type == MODEL_R2EFF:
807
808 if not has_exponential_exp_type():
809 raise RelaxError("The R2eff model with the fixed time period dispersion experiments cannot be optimised.")
810
811
812 minimise_r2eff(spins=spins, spin_ids=spin_ids, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iterations, constraints=constraints, scaling_matrix=scaling_matrix[model_index], verbosity=verbosity, sim_index=sim_index, lower=lower_i, upper=upper_i, inc=inc_i)
813
814
815 continue
816
817
818 command = Disp_minimise_command(spins=spins, spin_ids=spin_ids, sim_index=sim_index, scaling_matrix=scaling_matrix[model_index], min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iterations, constraints=constraints, verbosity=verbosity, lower=lower_i, upper=upper_i, inc=inc_i, fields=fields, param_names=get_param_names(spins=spins, full=True))
819
820
821 memo = Disp_memo(spins=spins, spin_ids=spin_ids, sim_index=sim_index, scaling_matrix=scaling_matrix[model_index], verbosity=verbosity)
822
823
824 processor.add_to_queue(command, memo)
825
826
828 """Return a description of the model.
829
830 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
831 @type model_info: list of SpinContainer instances, list of str
832 @return: The model description.
833 @rtype: str
834 """
835
836
837 return "The spin cluster %s." % model_info
838
839
841 """Loop over the spin groupings for one model applied to multiple spins.
842
843 @return: The list of spins per block will be yielded, as well as the list of spin IDs.
844 @rtype: tuple of list of SpinContainer instances and list of str
845 """
846
847
848 if not hasattr(cdp, 'model_type') or cdp.model_type == MODEL_R2EFF:
849
850 for spin, spin_id in spin_loop(return_id=True):
851
852 if not spin.select:
853 continue
854
855
856 yield [spin_id]
857
858
859 else:
860 for spin_ids in loop_cluster(skip_desel=False):
861 yield spin_ids
862
863
865 """Return the k, n, and chi2 model statistics.
866
867 k - number of parameters.
868 n - number of data points.
869 chi2 - the chi-squared value.
870
871
872 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
873 @type model_info: list of SpinContainer instances, list of str
874 @keyword spin_id: The spin ID string to override the model_info argument. This is ignored in the N-state model.
875 @type spin_id: None or str
876 @keyword global_stats: A parameter which determines if global or local statistics are returned. For the N-state model, this argument is ignored.
877 @type global_stats: None or bool
878 @return: The optimisation statistics, in tuple format, of the number of parameters (k), the number of data points (n), and the chi-squared value (chi2).
879 @rtype: tuple of (int, int, float)
880 """
881
882
883 if spin_id != None:
884 spins = [return_spin(spin_id)]
885 else:
886 spins = spin_ids_to_containers(model_info)
887
888
889 k = param_num(spins=spins)
890
891
892 n = 0
893 for spin in spins:
894
895 if not spin.select:
896 continue
897
898 n += len(spin.r2eff)
899
900
901 chi2 = None
902 for spin in spins:
903
904 if not spin.select:
905 continue
906
907 if hasattr(spin, 'chi2'):
908 chi2 = spin.chi2
909 break
910
911
912 return k, n, chi2
913
914
916 """Deselect spins which have insufficient data to support minimisation.
917
918 @keyword data_check: A flag to signal if the presence of base data is to be checked for.
919 @type data_check: bool
920 @keyword verbose: A flag which if True will allow printouts.
921 @type verbose: bool
922 """
923
924
925 check_mol_res_spin_data()
926 check_model_type()
927
928
929 proton_mmq_flag = has_proton_mmq_cpmg()
930
931
932 for spin, spin_id in spin_loop(return_id=True, skip_desel=True):
933
934 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H':
935 continue
936
937
938 proton = None
939 if proton_mmq_flag:
940
941 proton_spins = return_attached_protons(spin_id)
942
943
944 if len(proton_spins) > 1:
945 print("Multiple protons attached to the spin '%s', but one one attached proton is supported for the MMQ-type models." % spin_id)
946 spin.select = False
947 continue
948
949
950 if len(proton_spins):
951 proton = proton_spins[0]
952
953
954 if not hasattr(spin, 'r2eff') and not hasattr(proton, 'r2eff'):
955 print("No R2eff data could be found, deselecting the '%s' spin." % spin_id)
956 spin.select = False
957 continue
958
959
961 """Print out the model title.
962
963 @keyword prefix: The starting text of the title. This should be printed out first, followed by the model information text.
964 @type prefix: str
965 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
966 @type model_info: list of SpinContainer instances, list of str
967 """
968
969
970 subsection(file=sys.stdout, text=prefix+"the spin block %s"%model_info, prespace=2)
971
972
974 """Return the peak intensity data structure.
975
976 @param data_id: The spin ID string, as yielded by the base_data_loop() generator method.
977 @type data_id: str
978 @return: The peak intensity data structure.
979 @rtype: list of float
980 """
981
982
983 if cdp.model_type == MODEL_R2EFF:
984
985 spin, spin_id, exp_type, frq, offset, point = data_id
986
987
988 return spin.peak_intensity
989
990
991 else:
992 raise RelaxImplementError
993
994
996 """Return the standard deviation data structure.
997
998 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method.
999 @type data_id: SpinContainer instance and float
1000 @return: The standard deviation data structure.
1001 @rtype: list of float
1002 """
1003
1004
1005 if cdp.model_type == MODEL_R2EFF:
1006
1007 spin, spin_id, exp_type, frq, offset, point = data_id
1008
1009
1010 errors = []
1011 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point):
1012 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True))
1013
1014
1015 else:
1016
1017 spin, spin_id = data_id
1018
1019
1020 proton_mmq_flag = has_proton_mmq_cpmg()
1021
1022
1023 proton = None
1024 if proton_mmq_flag:
1025 proton = return_attached_protons(spin_id)[0]
1026
1027
1028 r2eff_err = {}
1029 if hasattr(spin, 'r2eff_err'):
1030 r2eff_err.update(spin.r2eff_err)
1031 if hasattr(proton, 'r2eff_err'):
1032 r2eff_err.update(proton.r2eff_err)
1033 return r2eff_err
1034
1035
1036 return errors
1037
1038
1040 """Return the standard deviation data structure, where standard deviation is from the overall gauss distribution described by the STD_fit of the goodness of fit, where STD_fit = sqrt(chi2/(N-p))
1041
1042 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method.
1043 @type data_id: SpinContainer instance and float
1044 @return: The standard deviation data structure.
1045 @rtype: list of float
1046 """
1047
1048
1049 if cdp.model_type == MODEL_R2EFF:
1050 raise RelaxError("Drawing errors from the gauss distribution described by the STD_fit of the goodness of fit, is not possible for the '%s' model."%MODEL_R2EFF)
1051
1052
1053 errors = self.return_error(data_id=data_id)
1054
1055
1056 spin, spin_id = data_id
1057
1058
1059 for spin_ids in self.model_loop():
1060
1061 if spin_id in spin_ids:
1062
1063 k, n, chi2 = self.model_statistics(model_info=spin_ids)
1064
1065
1066 dof = n - k
1067
1068
1069 red_chi2 = chi2 / float(dof)
1070
1071
1072 std_red_chi2 = sqrt(red_chi2)
1073
1074
1075 for id in errors:
1076 errors[id] = std_red_chi2
1077
1078 return errors
1079
1080
1082 """Return the value and error corresponding to the parameter.
1083
1084 If sim is set to an integer, return the value of the simulation and None.
1085
1086
1087 @param spin: The SpinContainer object.
1088 @type spin: SpinContainer
1089 @param param: The name of the parameter to return values for.
1090 @type param: str
1091 @keyword sim: The Monte Carlo simulation index.
1092 @type sim: None or int
1093 @keyword bc: The back-calculated data flag. If True, then the back-calculated data will be returned rather than the actual data.
1094 @type bc: bool
1095 @return: The value and error corresponding to
1096 @rtype: tuple of length 2 of floats or None
1097 """
1098
1099
1100 special_parameters = ['theta', 'w_eff']
1101
1102
1103 if param not in special_parameters:
1104 returnval = self._return_value_general(spin=spin, param=param, sim=sim, bc=bc)
1105 return returnval
1106
1107
1108
1109
1110 value = None
1111 error = None
1112
1113
1114 theta_spin_dic, Domega_spin_dic, w_eff_spin_dic, dic_key_list = calc_rotating_frame_params(spin=spin)
1115
1116
1117 if param == "theta":
1118 value = theta_spin_dic
1119
1120
1121 elif param == "w_eff":
1122 value = w_eff_spin_dic
1123
1124
1125 return value, error
1126
1127
1128 - def set_error(self, index, error, model_info=None):
1129 """Set the parameter errors.
1130
1131 @param index: The index of the parameter to set the errors for.
1132 @type index: int
1133 @param error: The error value.
1134 @type error: float
1135 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
1136 @type model_info: list of SpinContainer instances, list of str
1137 """
1138
1139
1140 spin_ids = model_info
1141 spins = spin_ids_to_containers(spin_ids)
1142
1143
1144 total_param_num = param_num(spins=spins)
1145
1146
1147 model_param = True
1148 if index >= total_param_num:
1149 model_param = False
1150
1151
1152 aux_params = []
1153 if 'pA' in spins[0].params:
1154 aux_params.append('pB')
1155 if 'pB' in spins[0].params:
1156 aux_params.append('pA')
1157 if 'kex' in spins[0].params:
1158 aux_params.append('tex')
1159 if 'tex' in spins[0].params:
1160 aux_params.append('kex')
1161
1162
1163 if model_param:
1164 param_name, si, mi = param_index_to_param_info(index=index, spins=spins)
1165 else:
1166 param_name = aux_params[index - total_param_num]
1167 si = None
1168 mi = None
1169
1170
1171 err_name = param_name + "_err"
1172
1173
1174 if param_name in ['r2eff', 'i0']:
1175
1176 if not hasattr(spins[si], err_name):
1177 setattr(spins[si], err_name, {})
1178
1179
1180 setattr(spins[si], err_name, error)
1181
1182
1183 else:
1184
1185 if si == None:
1186 for spin in spins:
1187 setattr(spin, err_name, error)
1188
1189
1190 else:
1191
1192 setattr(spins[si], err_name, error)
1193
1194
1195 - def set_param_values(self, param=None, value=None, index=None, spin_id=None, error=False, force=True):
1196 """Set the spin specific parameter values.
1197
1198 @keyword param: The parameter name list.
1199 @type param: list of str
1200 @keyword value: The parameter value list.
1201 @type value: list
1202 @keyword index: The index for parameters which are of the list-type.
1203 @type index: None or int
1204 @keyword spin_id: The spin identification string, only used for spin specific parameters.
1205 @type spin_id: None or str
1206 @keyword error: A flag which if True will allow the parameter errors to be set instead of the values.
1207 @type error: bool
1208 @keyword force: A flag which if True will cause current values to be overwritten. If False, a RelaxError will raised if the parameter value is already set.
1209 @type force: bool
1210 """
1211
1212
1213 is_str_list(param, 'parameter name')
1214 is_list(value, 'parameter value')
1215
1216
1217 model_index = -1
1218 for spin_ids in self.model_loop():
1219
1220 model_index += 1
1221
1222
1223 spins = spin_ids_to_containers(spin_ids)
1224
1225
1226 skip = True
1227 for spin in spins:
1228 if spin.select:
1229 skip = False
1230 if skip:
1231 continue
1232
1233
1234 for i in range(len(param)):
1235 param_i = param[i]
1236 value_i = value[i]
1237
1238
1239 if not self._PARAMS.contains(param_i):
1240 raise RelaxError("The parameter '%s' is not valid for this data pipe type." % param_i)
1241
1242
1243 if param_i in ['pA', 'pB', 'pC', 'kex', 'k_AB', 'k_BC', 'k_AC', 'tex', 'kB', 'kC', 'kex_AB', 'kex_BC', 'kex_AC']:
1244 selection_list = spin_ids
1245 else:
1246 selection_list = [spin_id]
1247
1248
1249 for selection in selection_list:
1250
1251 for spin in spin_loop(selection=selection):
1252
1253 obj_name = param_i
1254 if error:
1255 obj_name += '_err'
1256
1257
1258 if param_i in ['r1']:
1259
1260 for exp_type, frq, ei, mi in loop_exp_frq(return_indices=True):
1261
1262 key = generate_r20_key(exp_type=exp_type, frq=frq)
1263
1264
1265 if not hasattr(spin, obj_name):
1266 setattr(spin, obj_name, {})
1267
1268
1269 if index == None:
1270 obj = getattr(spin, obj_name)
1271 obj[key] = value_i
1272
1273
1274 elif mi == index:
1275 obj = getattr(spin, obj_name)
1276 obj[key] = value_i
1277
1278
1279 elif param_i in PARAMS_R20:
1280
1281 for exp_type, frq, ei, mi in loop_exp_frq(return_indices=True):
1282
1283 key = generate_r20_key(exp_type=exp_type, frq=frq)
1284
1285
1286 if not hasattr(spin, obj_name):
1287 setattr(spin, obj_name, {})
1288
1289
1290 if index == None:
1291 obj = getattr(spin, obj_name)
1292 obj[key] = value_i
1293
1294
1295 elif mi == index:
1296 obj = getattr(spin, obj_name)
1297 obj[key] = value_i
1298
1299
1300 elif param_i in ['r2eff', 'i0'] and not isinstance(value_i, dict):
1301
1302 for exp_type, frq, offset, point in loop_exp_frq_offset_point():
1303
1304 key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point)
1305
1306
1307 if not hasattr(spin, obj_name):
1308 setattr(spin, obj_name, {})
1309
1310
1311 obj = getattr(spin, obj_name)
1312 obj[key] = value_i
1313
1314
1315 else:
1316 setattr(spin, obj_name, value_i)
1317
1318
1320 """Set the simulation selection flag.
1321
1322 @param select_sim: The selection flag for the simulations.
1323 @type select_sim: bool
1324 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
1325 @type model_info: list of SpinContainer instances, list of str
1326 """
1327
1328
1329 spin_ids = model_info
1330 spins = spin_ids_to_containers(spin_ids)
1331
1332
1333 for spin in spins:
1334 spin.select_sim = deepcopy(select_sim)
1335
1336
1338 """Initialise the Monte Carlo parameter values."""
1339
1340
1341 param_names = self.data_names(set='params')
1342
1343
1344 pairs = {}
1345 pairs['kex'] = 'tex'
1346 pairs['tex'] = 'kex'
1347
1348
1349 pairs['pA'] = 'pB'
1350 pairs['pB'] = 'pA'
1351
1352
1353 pairs['k_AB'] = 'kex'
1354 pairs['k_BA'] = 'kex'
1355
1356
1357 min_names = self.data_names(set='min')
1358
1359
1360
1361
1362
1363
1364 for spin in spin_loop():
1365
1366 if not spin.select:
1367 continue
1368
1369
1370 for object_name in param_names:
1371
1372 sim_object_name = object_name + '_sim'
1373
1374
1375
1376
1377
1378
1379 for spin in spin_loop():
1380
1381 if not spin.select:
1382 continue
1383
1384
1385 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H':
1386 continue
1387
1388
1389 for object_name in param_names:
1390
1391 if not (object_name in spin.params or (object_name in pairs and pairs[object_name] in spin.params)):
1392 continue
1393
1394
1395 sim_object_name = object_name + '_sim'
1396
1397
1398 setattr(spin, sim_object_name, [])
1399
1400
1401 sim_object = getattr(spin, sim_object_name)
1402
1403
1404 for j in range(cdp.sim_number):
1405
1406 if object_name not in spin.params:
1407 value = deepcopy(getattr(spin, pairs[object_name]))
1408 else:
1409 value = deepcopy(getattr(spin, object_name))
1410
1411
1412 sim_object.append(value)
1413
1414
1415 for object_name in min_names:
1416
1417 sim_object_name = object_name + '_sim'
1418
1419
1420 setattr(spin, sim_object_name, [])
1421
1422
1423 sim_object = getattr(spin, sim_object_name)
1424
1425
1426 for j in range(cdp.sim_number):
1427
1428 sim_object.append(deepcopy(getattr(spin, object_name)))
1429
1430
1432 """Pack the Monte Carlo simulation data.
1433
1434 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method.
1435 @type data_id: SpinContainer instance and float
1436 @param sim_data: The Monte Carlo simulation data.
1437 @type sim_data: list of float
1438 """
1439
1440
1441 if cdp.model_type == MODEL_R2EFF:
1442
1443 spin, spin_id, exp_type, frq, offset, point = data_id
1444
1445
1446 if not hasattr(spin, 'peak_intensity_sim'):
1447 spin.peak_intensity_sim = {}
1448
1449
1450 ti = 0
1451 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point):
1452
1453 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)
1454
1455
1456 for int_key in int_keys:
1457
1458 if int_key in spin.peak_intensity_sim:
1459 raise RelaxError("Monte Carlo simulation data for the key '%s' already exists." % int_key)
1460
1461
1462 spin.peak_intensity_sim[int_key] = []
1463
1464
1465 for i in range(cdp.sim_number):
1466 spin.peak_intensity_sim[int_key].append(sim_data[i][ti])
1467
1468
1469 ti += 1
1470
1471
1472 else:
1473
1474 spin, spin_id = data_id
1475
1476
1477 proton_mmq_flag = has_proton_mmq_cpmg()
1478
1479
1480 proton = None
1481 if proton_mmq_flag:
1482 proton = return_attached_protons(spin_id)[0]
1483
1484
1485 spin.r2eff_sim = sim_data
1486 if proton != None:
1487 proton.r2eff_sim = sim_data
1488
1489
1491 """Return the array of simulation parameter values.
1492
1493 @param index: The index of the parameter to return the array of values for.
1494 @type index: int
1495 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
1496 @type model_info: list of SpinContainer instances, list of str
1497 @return: The array of simulation parameter values.
1498 @rtype: list of float
1499 """
1500
1501
1502 spin_ids = model_info
1503 spins = spin_ids_to_containers(spin_ids)
1504
1505
1506 total_param_num = param_num(spins=spins)
1507
1508
1509 model_param = True
1510 if index >= total_param_num:
1511 model_param = False
1512
1513
1514 aux_params = []
1515 for spin in spins:
1516 if not spin.select:
1517 continue
1518 if 'pA' in spin.params:
1519 aux_params.append('pB')
1520 if 'pB' in spin.params:
1521 aux_params.append('pA')
1522 if 'kex' in spin.params:
1523 aux_params.append('tex')
1524 if 'tex' in spin.params:
1525 aux_params.append('kex')
1526 break
1527
1528
1529 total_aux_num = total_param_num + len(aux_params)
1530 if index >= total_aux_num:
1531 return
1532
1533
1534 if model_param:
1535 param_name, si, mi = param_index_to_param_info(index=index, spins=spins)
1536 if not param_name in ['r2eff', 'i0']:
1537 if si == None:
1538 si = 0
1539
1540 else:
1541 param_name = aux_params[index - total_param_num]
1542 si = 0
1543 mi = 0
1544
1545
1546 sim_data = []
1547 if param_name == 'r2eff':
1548 for i in range(cdp.sim_number):
1549 sim_data.append(spins[si].r2eff_sim[i])
1550 elif param_name == 'i0':
1551 for i in range(cdp.sim_number):
1552 sim_data.append(spins[si].i0_sim[i])
1553
1554
1555 else:
1556 sim_data = getattr(spins[si], param_name + "_sim")
1557
1558
1559 if sim_data == []:
1560 sim_data = None
1561
1562
1563 return sim_data
1564
1565
1567 """Return the array of selected simulation flags.
1568
1569 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop().
1570 @type model_info: list of SpinContainer instances, list of str
1571 @return: The array of selected simulation flags.
1572 @rtype: list of int
1573 """
1574
1575
1576 spin_ids = model_info
1577 spins = spin_ids_to_containers(spin_ids)
1578
1579
1580 return spins[0].select_sim
1581