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