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