1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 """The Lipari-Szabo model-free analysis API object."""
25
26
27
28 import bmrblib
29 from copy import deepcopy
30 from math import pi
31 from minfx.grid import grid_split
32 from numpy import array, dot, float64, int32, zeros
33 from numpy.linalg import inv
34 from re import match, search
35 import string
36 import sys
37 from types import MethodType
38 from warnings import warn
39
40
41 from lib.arg_check import is_num_list, is_str_list
42 from lib.errors import RelaxError, RelaxFault, RelaxNoModelError, RelaxNoSequenceError, RelaxNoTensorError
43 from lib.float import isInf
44 from lib.periodic_table import periodic_table
45 from lib.physical_constants import h_bar, mu0
46 from lib.text.sectioning import subsection
47 from lib.warnings import RelaxDeselectWarning, RelaxWarning
48 from multi import Processor_box
49 from pipe_control import diffusion_tensor, interatomic, mol_res_spin, pipes, relax_data, sequence
50 from pipe_control.bmrb import list_sample_conditions
51 from pipe_control.exp_info import bmrb_write_citations, bmrb_write_methods, bmrb_write_software
52 from pipe_control.interatomic import return_interatom_list
53 from pipe_control.mol_res_spin import count_spins, exists_mol_res_spin_data, find_index, get_molecule_names, return_spin, return_spin_from_index, return_spin_indices, spin_loop
54 from pipe_control.spectrometer import check_spectrometer_setup
55 from specific_analyses.api_base import API_base
56 from specific_analyses.api_common import API_common
57 from specific_analyses.model_free.bmrb import sf_csa_read, sf_model_free_read, to_bmrb_model
58 from specific_analyses.model_free.data import compare_objects
59 from specific_analyses.model_free.molmol import Molmol
60 from specific_analyses.model_free.model import determine_model_type
61 from specific_analyses.model_free.parameters import are_mf_params_set, assemble_param_names, assemble_param_vector, linear_constraints
62 from specific_analyses.model_free.optimisation import MF_grid_command, MF_memo, MF_minimise_command, minimise_data_setup, relax_data_opt_structs
63 from specific_analyses.model_free.parameter_object import Model_free_params
64 from specific_analyses.model_free.pymol import Pymol
65 from target_functions.mf import Mf
66
67
69 """Parent class containing all the model-free specific functions."""
70
71
72 instance = None
73
93
94
95 - def back_calc_ri(self, spin_index=None, ri_id=None, ri_type=None, frq=None):
96 """Back-calculation of relaxation data from the model-free parameter values.
97
98 @keyword spin_index: The global spin index.
99 @type spin_index: int
100 @keyword ri_id: The relaxation data ID string.
101 @type ri_id: str
102 @keyword ri_type: The relaxation data type.
103 @type ri_type: str
104 @keyword frq: The field strength.
105 @type frq: float
106 @return: The back calculated relaxation data value corresponding to the index.
107 @rtype: float
108 """
109
110
111 spin, spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True)
112
113
114 if hasattr(cdp, 'diff_tensor') and (cdp.diff_tensor.type == 'spheroid' or cdp.diff_tensor.type == 'ellipsoid'):
115 interatoms = interatomic.return_interatom_list(spin_id)
116 for i in range(len(interatoms)):
117
118 if not interatoms[i].dipole_pair:
119 continue
120
121
122 if not hasattr(interatoms[i], 'vector') or interatoms[i].vector == None:
123 warn(RelaxDeselectWarning(spin_id, 'missing structural data'))
124 return
125
126
127 self.overfit_deselect(data_check=False, verbose=False)
128
129
130 value = self.minimise(min_algor='back_calc', min_options=(spin_index, ri_id, ri_type, frq))
131
132
133 return value
134
135
136 - def bmrb_read(self, file_path, version=None, sample_conditions=None):
174
175
177 """Write the model-free results to a BMRB NMR-STAR v3.1 formatted file.
178
179 @param file_path: The full file path.
180 @type file_path: str
181 @keyword version: The BMRB NMR-STAR dictionary format to output to.
182 @type version: str
183 """
184
185
186 check_spectrometer_setup(escalate=2)
187
188
189 cdp = pipes.get_pipe()
190
191
192 star = bmrblib.create_nmr_star('relax_model_free_results', file_path, version)
193
194
195 global_chi2 = None
196 if hasattr(cdp, 'chi2'):
197 global_chi2 = cdp.chi2
198
199
200 rex_frq = cdp.spectrometer_frq[cdp.ri_ids[0]]
201
202
203 mol_name_list = []
204 res_num_list = []
205 res_name_list = []
206 atom_name_list = []
207
208 csa_list = []
209 r_list = []
210 isotope_list = []
211 element_list = []
212
213 local_tm_list = []
214 s2_list = []
215 s2f_list = []
216 s2s_list = []
217 te_list = []
218 tf_list = []
219 ts_list = []
220 rex_list = []
221
222 local_tm_err_list = []
223 s2_err_list = []
224 s2f_err_list = []
225 s2s_err_list = []
226 te_err_list = []
227 tf_err_list = []
228 ts_err_list = []
229 rex_err_list = []
230
231 chi2_list = []
232 model_list = []
233
234
235 for spin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True):
236
237 if spin.name == 'H' or (hasattr(spin, 'element') and spin.element == 'H'):
238 warn(RelaxWarning("Skipping the proton spin '%s'." % spin_id))
239 continue
240
241
242 if res_num == None:
243 raise RelaxError("For the BMRB, the residue of spin '%s' must be numbered." % spin_id)
244 if res_name == None:
245 raise RelaxError("For the BMRB, the residue of spin '%s' must be named." % spin_id)
246 if spin.name == None:
247 raise RelaxError("For the BMRB, the spin '%s' must be named." % spin_id)
248 if not hasattr(spin, 'isotope') or spin.isotope == None:
249 raise RelaxError("For the BMRB, the spin isotope type of '%s' must be specified." % spin_id)
250 if not hasattr(spin, 'element') or spin.element == None:
251 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)
252
253
254 mol_name_list.append(mol_name)
255 res_num_list.append(res_num)
256 res_name_list.append(res_name)
257 atom_name_list.append(spin.name)
258
259
260 if hasattr(spin, 'csa'):
261 csa_list.append(spin.csa * 1e6)
262 else:
263 csa_list.append(None)
264
265
266 interatoms = return_interatom_list(spin_id)
267 for i in range(len(interatoms)):
268
269 if not interatoms[i].dipole_pair:
270 continue
271
272
273 if hasattr(interatoms[i], 'r'):
274 r_list.append(interatoms[i].r)
275 else:
276 r_list.append(None)
277
278
279 break
280
281
282 if hasattr(spin, 'isotope'):
283 isotope_list.append(int(spin.isotope.strip(string.ascii_letters)))
284 else:
285 isotope_list.append(None)
286
287
288 if hasattr(spin, 'element'):
289 element_list.append(spin.element)
290 else:
291 element_list.append(None)
292
293
294 if hasattr(spin, 'local_tm'):
295 local_tm_list.append(spin.local_tm)
296 else:
297 local_tm_list.append(None)
298 if hasattr(spin, 'local_tm_err'):
299 local_tm_err_list.append(spin.local_tm_err)
300 else:
301 local_tm_err_list.append(None)
302
303
304 s2_list.append(spin.s2)
305 s2f_list.append(spin.s2f)
306 s2s_list.append(spin.s2s)
307 te_list.append(spin.te)
308 tf_list.append(spin.tf)
309 ts_list.append(spin.ts)
310 if spin.rex == None:
311 rex_list.append(None)
312 else:
313 rex_list.append(spin.rex / (2.0*pi*rex_frq**2))
314
315
316 params = ['s2', 's2f', 's2s', 'te', 'tf', 'ts', 'rex']
317 for param in params:
318
319 err_list = locals()[param+'_err_list']
320
321
322 if hasattr(spin, param+'_err'):
323
324 val = getattr(spin, param+'_err')
325
326
327 if param == 'rex' and val != None:
328 val = val / (2.0*pi*rex_frq**2)
329
330
331 err_list.append(val)
332
333
334 else:
335 err_list.append(None)
336
337
338
339 chi2_list.append(spin.chi2)
340
341
342 model_list.append(to_bmrb_model(spin.model))
343
344
345 entity_ids = zeros(len(mol_name_list), int32)
346 mol_names = get_molecule_names()
347 for i in range(len(mol_name_list)):
348 for j in range(len(mol_names)):
349 if mol_name_list[i] == mol_names[j]:
350 entity_ids[i] = j+1
351
352
353
354
355
356
357 bmrb_write_citations(star)
358
359
360
361
362
363
364 mol_res_spin.bmrb_write_entity(star)
365
366
367
368
369
370
371 bmrb_write_methods(star)
372
373
374 software_ids, software_labels = bmrb_write_software(star)
375
376
377
378
379
380
381 star.chem_shift_anisotropy.add(entity_ids=entity_ids, res_nums=res_num_list, res_names=res_name_list, atom_names=atom_name_list, atom_types=element_list, isotope=isotope_list, csa=csa_list)
382
383
384
385
386
387
388 relax_data.bmrb_write(star)
389
390
391
392
393
394
395 for i in range(len(software_ids)):
396 if software_labels[i] == 'relax':
397 software_id = software_ids[i]
398
399
400 star.model_free.add(global_chi2=global_chi2, software_ids=[software_id], software_labels=['relax'], entity_ids=entity_ids, res_nums=res_num_list, res_names=res_name_list, atom_names=atom_name_list, atom_types=element_list, isotope=isotope_list, bond_length=r_list, local_tc=local_tm_list, s2=s2_list, s2f=s2f_list, s2s=s2s_list, te=te_list, tf=tf_list, ts=ts_list, rex=rex_list, local_tc_err=local_tm_err_list, s2_err=s2_err_list, s2f_err=s2f_err_list, s2s_err=s2s_err_list, te_err=te_err_list, tf_err=tf_err_list, ts_err=ts_err_list, rex_err=rex_err_list, rex_frq=rex_frq, chi2=chi2_list, model_fit=model_list)
401
402
403
404
405
406
407 diffusion_tensor.bmrb_write(star)
408
409
410
411 star.write()
412
413
414 - def calculate(self, spin_id=None, scaling_matrix=None, verbosity=1, sim_index=None):
415 """Calculation of the model-free chi-squared value.
416
417 @keyword spin_id: The spin identification string.
418 @type spin_id: str
419 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices.
420 @type scaling_matrix: list of numpy rank-2, float64 array or list of None
421 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
422 @type verbosity: int
423 @keyword sim_index: The optional MC simulation index.
424 @type sim_index: int
425 """
426
427
428 if not exists_mol_res_spin_data():
429 raise RelaxNoSequenceError
430
431
432 model_type = determine_model_type()
433
434
435 if model_type != 'local_tm' and not diffusion_tensor.diff_data_exists():
436 raise RelaxNoTensorError('diffusion')
437
438
439 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(cdp, 'structure'):
440 raise RelaxNoPdbError
441
442
443 for spin, id in spin_loop(spin_id, return_id=True):
444
445 if not spin.select:
446 continue
447
448
449 if not spin.model:
450 raise RelaxNoModelError
451
452
453 if not hasattr(spin, 'isotope'):
454 raise RelaxSpinTypeError
455
456
457 unset_param = are_mf_params_set(spin)
458 if unset_param != None:
459 raise RelaxNoValueError(unset_param)
460
461
462 if not hasattr(spin, 'csa') or spin.csa == None:
463 raise RelaxNoValueError("CSA")
464
465
466 interatoms = return_interatom_list(spin_id)
467 for interatom in interatoms:
468
469 if not interatom.dipole_pair:
470 continue
471
472
473 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(interatom, 'vector'):
474 raise RelaxNoVectorsError
475
476
477 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and hasattr(interatom, 'vector') and is_num_list(interatom.vector[0], raise_error=False):
478 raise RelaxMultiVectorError
479
480
481 if spin_id != interatom.spin_id1:
482 spin_id2 = interatom.spin_id1
483 else:
484 spin_id2 = interatom.spin_id2
485 spin2 = return_spin(spin_id2)
486
487
488 if not hasattr(spin2, 'isotope'):
489 raise RelaxSpinTypeError
490
491
492 if not hasattr(interatom, 'r') or interatom.r == None:
493 raise RelaxNoValueError("interatomic distance", spin_id=spin_id, spin_id2=spin_id2)
494
495
496 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'):
497 continue
498
499
500 for ri_id in cdp.ri_ids:
501
502 err = spin.ri_data_err[ri_id]
503
504
505 if err != None and err == 0.0:
506 raise RelaxError("Zero error for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (id, ri_id))
507 elif err != None and err < 0.0:
508 raise RelaxError("Negative error of %s for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (err, id, ri_id))
509
510
511 param_vector = assemble_param_vector(spin=spin, sim_index=sim_index)
512
513
514 data = relax_data_opt_structs(spin, sim_index=sim_index)
515
516
517 ri_data = [array(data[0])]
518 ri_data_err = [array(data[1])]
519 num_frq = [data[2]]
520 num_ri = [data[3]]
521 ri_labels = [data[4]]
522 frq = [data[5]]
523 remap_table = [data[6]]
524 noe_r1_table = [data[7]]
525 gx = [periodic_table.gyromagnetic_ratio(spin.isotope)]
526 if sim_index == None:
527 csa = [spin.csa]
528 else:
529 csa = [spin.csa_sim[sim_index]]
530
531
532 interatoms = return_interatom_list(spin_id)
533 for i in range(len(interatoms)):
534
535 if not interatoms[i].dipole_pair:
536 continue
537
538
539 if spin_id != interatoms[i].spin_id1:
540 spin_id2 = interatoms[i].spin_id1
541 else:
542 spin_id2 = interatoms[i].spin_id2
543 spin2 = return_spin(spin_id2)
544
545
546 if sim_index == None:
547 r = [interatoms[i].r]
548 else:
549 r = [interatoms[i].r_sim[sim_index]]
550
551
552 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere':
553 xh_unit_vectors = [interatoms[i].vector]
554 else:
555 xh_unit_vectors = [None]
556
557
558 gh = [periodic_table.gyromagnetic_ratio(spin2.isotope)]
559
560
561 num_params = [len(spin.params)]
562
563
564 param_values = [assemble_param_vector(model_type='mf')]
565
566
567 if model_type == 'local_tm':
568 diff_params = [spin.local_tm]
569 diff_type = 'sphere'
570 else:
571
572 diff_type = cdp.diff_tensor.type
573
574
575 if diff_type == 'sphere':
576 diff_params = [cdp.diff_tensor.tm]
577
578
579 elif diff_type == 'spheroid':
580 diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.theta, cdp.diff_tensor.phi]
581
582
583 elif diff_type == 'ellipsoid':
584 diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.Dr, cdp.diff_tensor.alpha, cdp.diff_tensor.beta, cdp.diff_tensor.gamma]
585
586
587 mf = Mf(init_params=param_vector, model_type='mf', diff_type=diff_type, diff_params=diff_params, num_spins=1, equations=[spin.equation], param_types=[spin.params], param_values=param_values, relax_data=ri_data, errors=ri_data_err, bond_length=r, csa=csa, num_frq=num_frq, frq=frq, num_ri=num_ri, remap_table=remap_table, noe_r1_table=noe_r1_table, ri_labels=ri_labels, gx=gx, gh=gh, h_bar=h_bar, mu0=mu0, num_params=num_params, vectors=xh_unit_vectors)
588
589
590 try:
591 chi2 = mf.func(param_vector)
592 except OverflowError:
593 chi2 = 1e200
594
595
596 if model_type == 'all' or model_type == 'diff':
597 cdp.chi2 = cdp.chi2 + chi2
598 else:
599 spin.chi2 = chi2
600
601
603 """Create the Monte Carlo Ri data.
604
605 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method.
606 @type data_id: str
607 @return: The Monte Carlo simulation data.
608 @rtype: list of floats
609 """
610
611
612 mc_data = []
613
614
615 spin = return_spin(data_id)
616 global_index = find_index(data_id)
617
618
619 if not spin.select:
620 return
621
622
623 if not hasattr(spin, 'model') or not spin.model:
624 return None
625
626
627 for ri_id in cdp.ri_ids:
628
629 value = self.back_calc_ri(spin_index=global_index, ri_id=ri_id, ri_type=cdp.ri_type[ri_id], frq=cdp.spectrometer_frq[ri_id])
630
631
632 mc_data.append(value)
633
634
635 return mc_data
636
637
639 """Initialise the spin specific data structures.
640
641 @param data: The spin ID string from the _base_data_loop_spin() method.
642 @type data: str
643 @keyword sim: The Monte Carlo simulation flag, which if true will initialise the simulation data structure.
644 @type sim: bool
645 """
646
647
648 spin = return_spin(data)
649
650
651 for name in self._PARAMS.loop(scope='spin'):
652
653 if name in ['ri_data', 'ri_data_bc', 'ri_data_err']:
654 continue
655
656
657 list_data = [ 'params' ]
658 if name in list_data:
659 init_data = []
660
661
662 init_data = None
663 if self._PARAMS.type(name) == bool:
664 init_data = False
665 if name == 'select':
666 init_data = True
667
668
669 if not hasattr(spin, name):
670 setattr(spin, name, init_data)
671
672
673 - def deselect(self, sim_index=None, model_info=None):
674 """Deselect models or simulations.
675
676 @keyword sim_index: The optional Monte Carlo simulation index. If None, then models will be deselected, otherwise the given simulation will.
677 @type sim_index: None or int
678 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
679 @type model_info: int
680 """
681
682
683 model_type = determine_model_type()
684
685
686 if model_type == 'mf' or model_type == 'local_tm':
687
688 spin = return_spin_from_index(model_info)
689
690
691 if sim_index == None:
692 spin.select = False
693
694
695 else:
696 spin.select_sim[sim_index] = False
697
698
699 else:
700
701 if sim_index == None:
702 raise RelaxError("Cannot deselect the global model.")
703
704
705 else:
706
707 cdp.select_sim[sim_index] = False
708
709
710 - def duplicate_data(self, pipe_from=None, pipe_to=None, model_info=None, global_stats=False, verbose=True):
711 """Duplicate the data specific to a single model-free model.
712
713 @keyword pipe_from: The data pipe to copy the data from.
714 @type pipe_from: str
715 @keyword pipe_to: The data pipe to copy the data to.
716 @type pipe_to: str
717 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
718 @type model_info: int
719 @keyword global_stats: The global statistics flag.
720 @type global_stats: bool
721 @keyword verbose: A flag which if True will cause info about each spin to be printed out as the sequence is generated.
722 @type verbose: bool
723 """
724
725
726 if model_info == None:
727 raise RelaxError("The model_info argument cannot be None.")
728
729
730 if not pipes.has_pipe(pipe_to):
731 pipes.create(pipe_to, pipe_type='mf', switch=False)
732
733
734 dp_from = pipes.get_pipe(pipe_from)
735 dp_to = pipes.get_pipe(pipe_to)
736
737
738 for data_name in dir(dp_from):
739
740 if data_name in ['diff_tensor', 'mol', 'interatomic', 'structure', 'exp_info']:
741 continue
742
743
744 if search('^_', data_name) or data_name in dp_from.__class__.__dict__:
745 continue
746
747
748 data_from = getattr(dp_from, data_name)
749
750
751 if hasattr(dp_to, data_name):
752
753 data_to = getattr(dp_to, data_name)
754
755
756 if data_from != data_to:
757 raise RelaxError("The object " + repr(data_name) + " is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
758
759
760 continue
761
762
763 setattr(dp_to, data_name, deepcopy(data_from))
764
765
766 if hasattr(dp_from, 'diff_tensor'):
767
768 if not hasattr(dp_to, 'diff_tensor'):
769 setattr(dp_to, 'diff_tensor', deepcopy(dp_from.diff_tensor))
770
771
772 else:
773
774 for data_name in dp_from.diff_tensor._mod_attr:
775
776 data_from = None
777 if hasattr(dp_from.diff_tensor, data_name):
778 data_from = getattr(dp_from.diff_tensor, data_name)
779
780
781 if data_from and not hasattr(dp_to.diff_tensor, data_name):
782 raise RelaxError("The diffusion tensor object " + repr(data_name) + " of the " + repr(pipe_from) + " data pipe is not located in the " + repr(pipe_to) + " data pipe.")
783 elif data_from:
784 data_to = getattr(dp_to.diff_tensor, data_name)
785 else:
786 continue
787
788
789 if data_from != data_to:
790 raise RelaxError("The object " + repr(data_name) + " is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
791
792
793 if hasattr(dp_from, 'structure'):
794
795 if not hasattr(dp_to, 'structure'):
796 setattr(dp_to, 'structure', deepcopy(dp_from.structure))
797
798
799 else:
800
801 compare_objects(dp_from.structure, dp_to.structure, pipe_from, pipe_to)
802
803
804 if len(dp_from.structure.structural_data) != len(dp_from.structure.structural_data):
805 raise RelaxError("The number of structural models is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
806
807
808 for i in range(len(dp_from.structure.structural_data)):
809
810 model_from = dp_from.structure.structural_data[i]
811 model_to = dp_to.structure.structural_data[i]
812
813
814 if model_from.num != model_to.num:
815 raise RelaxError("The structure models are not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
816
817
818 if len(model_from.mol) != len(model_to.mol):
819 raise RelaxError("The number of molecules is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".")
820
821
822 for mol_index in range(len(model_from.mol)):
823
824 compare_objects(model_from.mol[mol_index], model_to.mol[mol_index], pipe_from, pipe_to)
825
826
827 if dp_from.mol.is_empty():
828 return
829
830
831 if dp_to.mol.is_empty():
832 sequence.copy(pipe_from=pipe_from, pipe_to=pipe_to, preserve_select=True, verbose=verbose)
833
834
835 if dp_to.interatomic.is_empty():
836 interatomic.copy(pipe_from=pipe_from, pipe_to=pipe_to, verbose=verbose)
837
838
839 pipes.switch(pipe_from)
840 model_type = determine_model_type()
841
842
843 spin, spin_id = return_spin_from_index(model_info, pipe=pipe_from, return_spin_id=True)
844 if model_type == 'mf' or (model_type == 'local_tm' and not global_stats):
845
846 for name in dir(spin):
847
848 if search('^__', name):
849 continue
850
851
852 obj = getattr(spin, name)
853
854
855 if isinstance(obj, MethodType):
856 continue
857
858
859 new_obj = deepcopy(getattr(spin, name))
860 setattr(dp_to.mol[spin._mol_index].res[spin._res_index].spin[spin._spin_index], name, new_obj)
861
862
863 interatoms = interatomic.return_interatom_list(spin_id)
864 for interatom in interatoms:
865
866 if not interatom.dipole_pair:
867 continue
868
869
870 if spin_id != interatom.spin_id1:
871 spin_id2 = interatom.spin_id1
872 else:
873 spin_id2 = interatom.spin_id2
874 spin2 = return_spin(spin_id2)
875
876
877 mol_index, res_index, spin_index = return_spin_indices(spin_id2)
878 dp_to.mol[mol_index].res[res_index].spin[spin_index] = deepcopy(spin2)
879
880
881 else:
882
883 dp_to.mol = deepcopy(dp_from.mol)
884
885
886 - def eliminate(self, name, value, args, sim=None, model_info=None):
887 """Model-free model elimination, parameter by parameter.
888
889 @param name: The parameter name.
890 @type name: str
891 @param value: The parameter value.
892 @type value: float
893 @param args: The c1 and c2 elimination constant overrides.
894 @type args: None or tuple of float
895 @keyword sim: The Monte Carlo simulation index.
896 @type sim: int
897 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
898 @type model_info: int
899 @return: True if the model is to be eliminated, False otherwise.
900 @rtype: bool
901 """
902
903
904 c1 = 50.0 * 1e-9
905 c2 = 1.5
906
907
908 if args != None:
909 c1, c2 = args
910
911
912 model_type = determine_model_type()
913
914
915 if model_type != 'mf' and model_type != 'local_tm':
916 raise RelaxError("Elimination of the global model is not yet supported.")
917
918
919 spin, spin_id = return_spin_from_index(model_info, return_spin_id=True)
920
921
922 if model_type == 'local_tm':
923 tm = spin.local_tm
924 else:
925 tm = cdp.diff_tensor.tm
926
927
928 if tm == None:
929 return False
930
931
932 if name == 'local_tm' and value >= c1:
933 if sim == None:
934 print("Data pipe '%s': The local tm parameter of %.5g is greater than %.5g, eliminating spin system '%s'." % (pipes.cdp_name(), value, c1, spin_id))
935 else:
936 print("Data pipe '%s': The local tm parameter of %.5g is greater than %.5g, eliminating simulation %i of spin system '%s'." % (pipes.cdp_name(), value, c1, sim, spin_id))
937 return True
938
939
940 if match('t[efs]', name) and value >= c2 * tm:
941 if sim == None:
942 print("Data pipe '%s': The %s value of %.5g is greater than %.5g, eliminating spin system '%s'." % (pipes.cdp_name(), name, value, c2*tm, spin_id))
943 else:
944 print("Data pipe '%s': The %s value of %.5g is greater than %.5g, eliminating simulation %i of spin system '%s'." % (pipes.cdp_name(), name, value, c2*tm, sim, spin_id))
945 return True
946
947
948 return False
949
950
972
973
1006
1007
1008 - def grid_search(self, lower=None, upper=None, inc=None, scaling_matrix=None, constraints=True, verbosity=1, sim_index=None):
1009 """The model-free grid search function.
1010
1011 @keyword lower: The per-model lower bounds of the grid search which must be equal to the number of parameters in the model.
1012 @type lower: list of lists of numbers
1013 @keyword upper: The per-model upper bounds of the grid search which must be equal to the number of parameters in the model.
1014 @type upper: list of lists of numbers
1015 @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.
1016 @type inc: list of lists of int
1017 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices.
1018 @type scaling_matrix: list of numpy rank-2, float64 array or list of None
1019 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used.
1020 @type constraints: bool
1021 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity.
1022 @type verbosity: int
1023 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised.
1024 @type sim_index: int
1025 """
1026
1027
1028 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
1029
1030
1032 """Create bounds for the OpenDX mapping function.
1033
1034 @param param: The name of the parameter to return the lower and upper bounds of.
1035 @type param: str
1036 @param spin_id: The spin identification string.
1037 @type spin_id: str
1038 @return: The upper and lower bounds of the parameter.
1039 @rtype: list of float
1040 """
1041
1042
1043 if self._PARAMS.scope(param) == 'global':
1044 return diffusion_tensor.map_bounds(param)
1045
1046
1047 spin = return_spin(spin_id)
1048
1049
1050 if search('^s2', param):
1051 return [0.0, 1.0]
1052
1053
1054 elif search('^t', param) or param == 'local_tm':
1055 return [0.0, 1e-8]
1056
1057
1058 elif param == 'rex':
1059 return [0.0, 30.0 / (2.0 * pi * cdp.spectrometer_frq[cdp.ri_ids[0]])**2]
1060
1061
1062 elif param == 'r':
1063 return [1.0 * 1e-10, 1.1 * 1e-10]
1064
1065
1066 elif param == 'csa':
1067 return [-100 * 1e-6, -300 * 1e-6]
1068
1069
1070 - 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):
1071 """Model-free minimisation function.
1072
1073 Three categories of models exist for which the approach to minimisation is different. These
1074 are:
1075
1076 Single spin optimisations: The 'mf' and 'local_tm' model types which are the
1077 model-free parameters for single spins, optionally with a local tm parameter. These
1078 models have no global parameters.
1079
1080 Diffusion tensor optimisations: The 'diff' diffusion tensor model type. No spin
1081 specific parameters exist.
1082
1083 Optimisation of everything: The 'all' model type consisting of all model-free and all
1084 diffusion tensor parameters.
1085
1086
1087 @keyword min_algor: The minimisation algorithm to use.
1088 @type min_algor: str
1089 @keyword min_options: An array of options to be used by the minimisation algorithm.
1090 @type min_options: array of str
1091 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
1092 @type func_tol: None or float
1093 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check.
1094 @type grad_tol: None or float
1095 @keyword max_iterations: The maximum number of iterations for the algorithm.
1096 @type max_iterations: int
1097 @keyword constraints: If True, constraints are used during optimisation.
1098 @type constraints: bool
1099 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices.
1100 @type scaling_matrix: list of numpy rank-2, float64 array or list of None
1101 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity.
1102 @type verbosity: int
1103 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired.
1104 @type sim_index: None or int
1105 @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.
1106 @type lower: list of lists of numbers
1107 @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.
1108 @type upper: list of lists of numbers
1109 @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.
1110 @type inc: list of lists of int
1111 """
1112
1113
1114 if not exists_mol_res_spin_data():
1115 raise RelaxNoSequenceError
1116
1117
1118 for spin in spin_loop():
1119
1120 if not spin.select:
1121 continue
1122
1123
1124 if not spin.model:
1125 raise RelaxNoModelError
1126
1127
1128 if not hasattr(spin, 'isotope'):
1129 raise RelaxSpinTypeError
1130
1131
1132 data_store = Data_container()
1133 opt_params = Data_container()
1134
1135
1136 data_store.h_bar = h_bar
1137 data_store.mu0 = mu0
1138 opt_params.min_algor = min_algor
1139 opt_params.min_options = min_options
1140 opt_params.func_tol = func_tol
1141 opt_params.grad_tol = grad_tol
1142 opt_params.max_iterations = max_iterations
1143
1144
1145 opt_params.verbosity = verbosity
1146
1147
1148 data_store.model_type = determine_model_type()
1149 if not data_store.model_type:
1150 return
1151
1152
1153 if min_algor == 'back_calc' and data_store.model_type != 'local_tm':
1154 data_store.model_type = 'mf'
1155
1156
1157 if data_store.model_type != 'local_tm' and not diffusion_tensor.diff_data_exists():
1158 raise RelaxNoTensorError('diffusion')
1159
1160
1161 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere':
1162
1163 if not hasattr(cdp, 'structure'):
1164 raise RelaxNoPdbError
1165
1166
1167 for spin, spin_id in spin_loop(return_id=True):
1168
1169 if not spin.select:
1170 continue
1171
1172
1173 interatoms = return_interatom_list(spin_id)
1174
1175
1176 for i in range(len(interatoms)):
1177
1178 if not interatoms[i].dipole_pair:
1179 continue
1180
1181
1182 if not hasattr(interatoms[i], 'vector'):
1183 raise RelaxNoVectorsError
1184
1185
1186 if data_store.model_type == 'diff':
1187
1188 for spin in spin_loop():
1189 unset_param = are_mf_params_set(spin)
1190 if unset_param != None:
1191 raise RelaxNoValueError(unset_param)
1192
1193
1194 if verbosity >= 1:
1195 if data_store.model_type == 'mf':
1196 print("Only the model-free parameters for single spins will be used.")
1197 elif data_store.model_type == 'local_mf':
1198 print("Only a local tm value together with the model-free parameters for single spins will be used.")
1199 elif data_store.model_type == 'diff':
1200 print("Only diffusion tensor parameters will be used.")
1201 elif data_store.model_type == 'all':
1202 print("The diffusion tensor parameters together with the model-free parameters for all spins will be used.")
1203
1204
1205 for spin, spin_id in spin_loop(return_id=True):
1206
1207 if not spin.select:
1208 continue
1209
1210
1211 if not hasattr(spin, 'csa') or spin.csa == None:
1212 raise RelaxNoValueError("CSA")
1213
1214
1215 interatoms = return_interatom_list(spin_id)
1216
1217
1218 count = 0
1219 for i in range(len(interatoms)):
1220
1221 if not interatoms[i].dipole_pair:
1222 continue
1223
1224
1225 if not hasattr(interatoms[i], 'r') or interatoms[i].r == None:
1226 raise RelaxNoValueError("interatomic distance", spin_id=interatoms[i].spin_id1, spin_id2=interatoms[i].spin_id2)
1227
1228
1229 count += 1
1230
1231
1232 if count > 1:
1233 raise RelaxError("The spin '%s' has %s dipolar relaxation interactions defined, but only a maximum of one is currently supported." % (spin_id, count))
1234
1235
1236 if data_store.model_type == 'mf' or data_store.model_type == 'local_tm':
1237 num_data_sets = 1
1238 data_store.num_spins = 1
1239 elif data_store.model_type == 'diff' or data_store.model_type == 'all':
1240 num_data_sets = count_spins(skip_desel=False)
1241 data_store.num_spins = count_spins()
1242
1243
1244 if min_algor == 'back_calc':
1245 num_data_sets = 0
1246 data_store.num_spins = 1
1247
1248
1249 processor_box = Processor_box()
1250 processor = processor_box.processor
1251
1252
1253 for index in self.model_loop():
1254
1255 if data_store.model_type == 'diff' or data_store.model_type == 'all':
1256 spin_index = None
1257 spin, data_store.spin_id = None, None
1258 elif min_algor == 'back_calc':
1259 spin_index = opt_params.min_options[0]
1260 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True)
1261 else:
1262 spin_index = index
1263 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True)
1264
1265
1266 if spin and (data_store.model_type == 'mf' or data_store.model_type == 'local_tm') and not min_algor == 'back_calc':
1267
1268 if not spin.select:
1269 continue
1270
1271
1272 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'):
1273 continue
1274
1275
1276 if spin and (data_store.model_type == 'mf' or data_store.model_type == 'local_tm'):
1277 interatoms = return_interatom_list(data_store.spin_id)
1278 if not len(interatoms):
1279 continue
1280
1281
1282 if min_algor == 'back_calc':
1283
1284 opt_params.param_vector = assemble_param_vector(spin=spin, model_type=data_store.model_type)
1285
1286
1287 data_store.scaling_matrix = None
1288
1289 else:
1290
1291 opt_params.param_vector = assemble_param_vector(spin=spin, sim_index=sim_index)
1292
1293
1294 num_params = len(opt_params.param_vector)
1295
1296
1297 data_store.scaling_matrix = scaling_matrix[index]
1298 if data_store.scaling_matrix != None:
1299 opt_params.param_vector = dot(inv(data_store.scaling_matrix), opt_params.param_vector)
1300
1301
1302 opt_params.lower, opt_params.upper, opt_params.inc = None, None, None
1303 if lower != None:
1304 opt_params.lower = lower[index]
1305 if upper != None:
1306 opt_params.upper = upper[index]
1307 if inc != None:
1308 opt_params.inc = inc[index]
1309
1310
1311 if match('^[Ss]et', min_algor):
1312 opt_params.min_options = dot(inv(data_store.scaling_matrix), opt_params.min_options)
1313
1314
1315 if constraints:
1316 opt_params.A, opt_params.b = linear_constraints(num_params, model_type=data_store.model_type, spin=spin, scaling_matrix=data_store.scaling_matrix)
1317 else:
1318 opt_params.A, opt_params.b = None, None
1319
1320
1321 minimise_data_setup(data_store, min_algor, num_data_sets, opt_params.min_options, spin=spin, sim_index=sim_index)
1322
1323
1324 if constraints and not match('^[Gg]rid', min_algor):
1325 algor = opt_params.min_options[0]
1326 else:
1327 algor = min_algor
1328
1329
1330 if min_algor == 'back_calc' or match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor):
1331 mf = Mf(init_params=opt_params.param_vector, model_type=data_store.model_type, diff_type=data_store.diff_type, diff_params=data_store.diff_params, scaling_matrix=data_store.scaling_matrix, num_spins=data_store.num_spins, equations=data_store.equations, param_types=data_store.param_types, param_values=data_store.param_values, relax_data=data_store.ri_data, errors=data_store.ri_data_err, bond_length=data_store.r, csa=data_store.csa, num_frq=data_store.num_frq, frq=data_store.frq, num_ri=data_store.num_ri, remap_table=data_store.remap_table, noe_r1_table=data_store.noe_r1_table, ri_labels=data_store.ri_types, gx=data_store.gx, gh=data_store.gh, h_bar=data_store.h_bar, mu0=data_store.mu0, num_params=data_store.num_params, vectors=data_store.xh_unit_vectors)
1332
1333
1334 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor):
1335
1336 number_ri = 0
1337 for k in range(len(ri_data_err)):
1338 number_ri = number_ri + len(ri_data_err[k])
1339
1340
1341 lm_error = zeros(number_ri, float64)
1342 index = 0
1343 for k in range(len(ri_data_err)):
1344 lm_error[index:index+len(ri_data_err[k])] = ri_data_err[k]
1345 index = index + len(ri_data_err[k])
1346
1347 opt_params.min_options = opt_params.min_options + (mf.lm_dri, lm_error)
1348
1349
1350 if min_algor == 'back_calc':
1351 return mf.calc_ri()
1352
1353
1354 if match('^[Gg]rid', min_algor) and data_store.model_type == 'diff':
1355
1356 print("Parallelised diffusion tensor grid search.")
1357
1358
1359 for subdivision in grid_split(divisions=processor.processor_size(), lower=opt_params.lower, upper=opt_params.upper, inc=opt_params.inc, verbosity=verbosity):
1360
1361 opt_params.subdivision = subdivision
1362
1363
1364 command = MF_grid_command()
1365
1366
1367 command.store_data(deepcopy(data_store), deepcopy(opt_params))
1368
1369
1370 memo = MF_memo(model_free=self, model_type=data_store.model_type, spin=spin, sim_index=sim_index, scaling_matrix=data_store.scaling_matrix)
1371 processor.add_to_queue(command, memo)
1372
1373
1374 processor.run_queue()
1375
1376
1377 return
1378
1379
1380 if search('^[Gg]rid', min_algor):
1381 command = MF_grid_command()
1382
1383
1384 else:
1385 command = MF_minimise_command()
1386
1387
1388 command.store_data(deepcopy(data_store), deepcopy(opt_params))
1389
1390
1391 memo = MF_memo(model_free=self, model_type=data_store.model_type, spin=spin, sim_index=sim_index, scaling_matrix=data_store.scaling_matrix)
1392 processor.add_to_queue(command, memo)
1393
1394
1395 processor.run_queue()
1396
1397
1399 """Return a description of the model.
1400
1401 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
1402 @type model_info: int
1403 @return: The model description.
1404 @rtype: str
1405 """
1406
1407
1408 model_type = determine_model_type()
1409
1410
1411 if model_type == 'all':
1412 return "Global model - all diffusion tensor parameters and spin specific model-free parameters."
1413 elif model_type == 'diff':
1414 return "Diffusion tensor model."
1415
1416
1417 else:
1418
1419 spin, spin_id = return_spin_from_index(model_info, return_spin_id=True)
1420
1421
1422 return "Model-free model of spin '%s'." % spin_id
1423
1424
1426 """Generator method for looping over the models (global or local).
1427
1428 If the model type is 'all' or 'diff', then this yields the single value of zero. Otherwise
1429 the global spin index is yielded.
1430
1431
1432 @return: The model index. This index is zero for the global models or equal to the global spin
1433 index (which covers the molecule, residue, and spin indices).
1434 @rtype: int
1435 """
1436
1437
1438 model_type = determine_model_type()
1439
1440
1441 if model_type == 'all' or model_type == 'diff':
1442 yield 0
1443
1444
1445 else:
1446
1447 global_index = -1
1448 for spin in spin_loop():
1449
1450 global_index = global_index + 1
1451
1452
1453 yield global_index
1454
1455
1457 """Return the k, n, and chi2 model statistics.
1458
1459 k - number of parameters.
1460 n - number of data points.
1461 chi2 - the chi-squared value.
1462
1463
1464 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
1465 @type model_info: int
1466 @keyword spin_id: The spin identification string. Either this or the instance keyword argument must be supplied.
1467 @type spin_id: None or str
1468 @keyword global_stats: A parameter which determines if global or local statistics are returned. If None, then the appropriateness of global or local statistics is automatically determined.
1469 @type global_stats: None or bool
1470 @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).
1471 @rtype: tuple of (int, int, float)
1472 """
1473
1474
1475 if model_info == None and spin_id == None:
1476 raise RelaxError("Either the model_info or spin_id argument must be supplied.")
1477 elif model_info != None and spin_id != None:
1478 raise RelaxError("The model_info arg " + repr(model_info) + " and spin_id arg " + repr(spin_id) + " clash. Only one should be supplied.")
1479
1480
1481 model_type = determine_model_type()
1482
1483
1484 if global_stats == None:
1485 if model_type in ['mf', 'local_tm']:
1486 global_stats = False
1487 else:
1488 global_stats = True
1489
1490
1491 if not global_stats:
1492
1493 if spin_id:
1494 spin = return_spin(spin_id)
1495 else:
1496 spin = return_spin_from_index(model_info)
1497
1498
1499 if not spin.select:
1500 return None, None, None
1501
1502
1503 if not hasattr(spin, 'ri_data'):
1504 return None, None, None
1505
1506
1507 param_vector = assemble_param_vector(spin=spin)
1508 k = len(param_vector)
1509
1510
1511 n = len(spin.ri_data)
1512
1513
1514 chi2 = spin.chi2
1515
1516
1517 elif global_stats:
1518
1519 param_vector = assemble_param_vector()
1520 k = len(param_vector)
1521
1522
1523 n = 0
1524 chi2 = 0
1525 for spin in spin_loop():
1526
1527 if not spin.select:
1528 continue
1529
1530
1531 if not hasattr(spin, 'ri_data') or not len(spin.ri_data):
1532 continue
1533
1534 n = n + len(spin.ri_data)
1535
1536
1537 if model_type == 'local_tm':
1538 chi2 = chi2 + spin.chi2
1539
1540
1541 if model_type != 'local_tm':
1542 if not hasattr(cdp, 'chi2'):
1543 raise RelaxError("Global statistics are not available, most likely because the global model has not been optimised.")
1544 chi2 = cdp.chi2
1545
1546
1547 return k, n, chi2
1548
1549
1551 """Return the type of the model, either being 'local' or 'global'.
1552
1553 @return: The model type, one of 'local' or 'global'.
1554 @rtype: str
1555 """
1556
1557
1558 model_type = determine_model_type()
1559
1560
1561 if model_type in ['all', 'diff']:
1562 return 'global'
1563
1564
1565 else:
1566 return 'local'
1567
1568
1570 """Function for returning the number of instances.
1571
1572 @return: The number of instances used for optimisation. Either the number of spins if
1573 the local optimisations are setup ('mf' and 'local_tm'), or 1 for the global
1574 models.
1575 @rtype: int
1576 """
1577
1578
1579 if not exists_mol_res_spin_data():
1580 return 0
1581
1582
1583 model_type = determine_model_type()
1584
1585
1586 if model_type == 'mf' or model_type == 'local_tm':
1587 return count_spins()
1588
1589
1590 elif model_type == 'diff' or model_type == 'all':
1591 return 1
1592
1593
1594 else:
1595 raise RelaxFault
1596
1597
1599 """Deselect spins which have insufficient data to support minimisation.
1600
1601 @keyword data_check: A flag to signal if the presence of base data is to be checked for.
1602 @type data_check: bool
1603 @keyword verbose: A flag which if True will allow printouts.
1604 @type verbose: bool
1605 """
1606
1607
1608 if verbose:
1609 print("\nOver-fit spin deselection:")
1610
1611
1612 if not exists_mol_res_spin_data():
1613 raise RelaxNoSequenceError
1614
1615
1616 need_vect = False
1617 if hasattr(cdp, 'diff_tensor') and (cdp.diff_tensor.type == 'spheroid' or cdp.diff_tensor.type == 'ellipsoid'):
1618 need_vect = True
1619
1620
1621 deselect_flag = False
1622 spin_count = 0
1623 for spin, spin_id in spin_loop(return_id=True):
1624
1625 if not spin.select:
1626 continue
1627
1628
1629 interatoms = interatomic.return_interatom_list(spin_id)
1630
1631
1632 dipole_relax = False
1633 for i in range(len(interatoms)):
1634
1635 if not interatoms[i].dipole_pair:
1636 continue
1637
1638
1639 if spin_id != interatoms[i].spin_id1:
1640 spin_id2 = interatoms[i].spin_id1
1641 else:
1642 spin_id2 = interatoms[i].spin_id2
1643 spin2 = return_spin(spin_id2)
1644
1645
1646 dipole_relax = True
1647
1648
1649 if not dipole_relax or not hasattr(spin, 'csa') or spin.csa == None:
1650 warn(RelaxDeselectWarning(spin_id, 'an absence of relaxation mechanisms'))
1651 spin.select = False
1652 deselect_flag = True
1653 continue
1654
1655
1656 if data_check:
1657
1658 data_points = 0
1659 inf_data = False
1660 if hasattr(cdp, 'ri_ids') and hasattr(spin, 'ri_data'):
1661 for id in cdp.ri_ids:
1662 if id in spin.ri_data and spin.ri_data[id] != None:
1663 data_points += 1
1664
1665
1666 if isInf(spin.ri_data[id]):
1667 inf_data = True
1668
1669
1670 if inf_data:
1671 warn(RelaxDeselectWarning(spin_id, 'infinite relaxation data'))
1672 spin.select = False
1673 deselect_flag = True
1674 continue
1675
1676
1677 if not hasattr(spin, 'ri_data'):
1678 warn(RelaxDeselectWarning(spin_id, 'missing relaxation data'))
1679 spin.select = False
1680 deselect_flag = True
1681 continue
1682
1683
1684 elif data_points < 3:
1685 warn(RelaxDeselectWarning(spin_id, 'insufficient relaxation data, 3 or more data points are required'))
1686 spin.select = False
1687 deselect_flag = True
1688 continue
1689
1690
1691 elif hasattr(spin, 'params') and spin.params and len(spin.params) > data_points:
1692 warn(RelaxDeselectWarning(spin_id, 'over-fitting - more parameters than data points'))
1693 spin.select = False
1694 deselect_flag = True
1695 continue
1696
1697
1698 for i in range(len(interatoms)):
1699
1700 if not interatoms[i].dipole_pair:
1701 continue
1702
1703
1704 if need_vect:
1705 if not hasattr(interatoms[i], 'vector') or interatoms[i].vector == None:
1706 warn(RelaxDeselectWarning(spin_id, 'missing structural data'))
1707 spin.select = False
1708 deselect_flag = True
1709 continue
1710
1711
1712 spin_count += 1
1713
1714
1715 if spin_count == 0:
1716 warn(RelaxWarning("No spins are selected therefore the optimisation or calculation cannot proceed."))
1717
1718
1719 if verbose and not deselect_flag:
1720 print("No spins have been deselected.")
1721
1722
1724 """Print out the model title.
1725
1726 @keyword prefix: The starting text of the title. This should be printed out first, followed by the model information text.
1727 @type prefix: str
1728 @keyword model_info: The model information from model_loop().
1729 @type model_info: unknown
1730 """
1731
1732
1733 model_type = determine_model_type()
1734
1735
1736 if model_type == 'mf' or model_type == 'local_tm':
1737 spin, spin_id = return_spin_from_index(global_index=model_info, return_spin_id=True)
1738 text = "%sSpin '%s'" % (prefix, spin_id)
1739
1740
1741 else:
1742 text = prefix + "Global model"
1743
1744
1745 subsection(file=sys.stdout, text=text, prespace=2)
1746
1747
1748 - def set_error(self, index, error, model_info=None):
1749 """Set the parameter errors.
1750
1751 @param index: The index of the parameter to set the errors for.
1752 @type index: int
1753 @param error: The error value.
1754 @type error: float
1755 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
1756 @type model_info: int
1757 """
1758
1759
1760 inc = 0
1761
1762
1763 model_type = determine_model_type()
1764
1765
1766 param_names = self.data_names(set='params', scope='spin')
1767
1768
1769
1770
1771
1772 if model_type == 'diff' or model_type == 'all':
1773
1774 if cdp.diff_tensor.type == 'sphere':
1775
1776 if index == 0:
1777 cdp.diff_tensor.set(param='tm', value=error, category='err')
1778
1779
1780 inc = inc + 1
1781
1782
1783 elif cdp.diff_tensor.type == 'spheroid':
1784
1785 if index == 0:
1786 cdp.diff_tensor.set(param='tm', value=error, category='err')
1787 elif index == 1:
1788 cdp.diff_tensor.set(param='Da', value=error, category='err')
1789 elif index == 2:
1790 cdp.diff_tensor.set(param='theta', value=error, category='err')
1791 elif index == 3:
1792 cdp.diff_tensor.set(param='phi', value=error, category='err')
1793
1794
1795 inc = inc + 4
1796
1797
1798 elif cdp.diff_tensor.type == 'ellipsoid':
1799
1800 if index == 0:
1801 cdp.diff_tensor.set(param='tm', value=error, category='err')
1802 elif index == 1:
1803 cdp.diff_tensor.set(param='Da', value=error, category='err')
1804 elif index == 2:
1805 cdp.diff_tensor.set(param='Dr', value=error, category='err')
1806 elif index == 3:
1807 cdp.diff_tensor.set(param='alpha', value=error, category='err')
1808 elif index == 4:
1809 cdp.diff_tensor.set(param='beta', value=error, category='err')
1810 elif index == 5:
1811 cdp.diff_tensor.set(param='gamma', value=error, category='err')
1812
1813
1814 inc = inc + 6
1815
1816
1817
1818
1819
1820 if model_type == 'all':
1821
1822 for spin in spin_loop():
1823
1824 if not spin.select:
1825 continue
1826
1827
1828 for param in param_names:
1829
1830 if index == inc:
1831 setattr(spin, param + "_err", error)
1832
1833
1834 inc = inc + 1
1835
1836
1837
1838
1839
1840 if model_type == 'mf' or model_type == 'local_tm':
1841
1842 spin = return_spin_from_index(model_info)
1843
1844
1845 if not spin.select:
1846 return
1847
1848
1849 for param in param_names:
1850
1851 if index == inc:
1852 setattr(spin, param + "_err", error)
1853
1854
1855 inc = inc + 1
1856
1857
1858 - def set_param_values(self, param=None, value=None, index=None, spin_id=None, error=False, force=True):
1859 """Set the model-free parameter values.
1860
1861 @keyword param: The parameter name list.
1862 @type param: list of str
1863 @keyword value: The parameter value list.
1864 @type value: list
1865 @keyword index: The index for parameters which are of the list-type. This is unused.
1866 @type index: None or int
1867 @keyword spin_id: The spin identification string, only used for spin specific parameters.
1868 @type spin_id: None or str
1869 @keyword error: A flag which if True will allow the parameter errors to be set instead of the values.
1870 @type error: bool
1871 @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.
1872 @type force: bool
1873 """
1874
1875
1876 is_str_list(param, 'parameter name')
1877
1878
1879 diff_params = []
1880 diff_vals = []
1881 mf_params = []
1882 mf_vals = []
1883 for i in range(len(param)):
1884
1885 if self._PARAMS.scope(param[i]) == 'global':
1886 if error:
1887 diff_params.append(param[i] + '_err')
1888 else:
1889 diff_params.append(param[i])
1890 diff_vals.append(value[i])
1891
1892
1893 else:
1894 mf_params.append(param[i])
1895 mf_vals.append(value[i])
1896
1897
1898 if diff_params:
1899 diffusion_tensor.set(value=diff_vals, param=diff_params)
1900
1901
1902 for i in range(len(mf_params)):
1903
1904 if mf_params[i] not in self.data_names(set='params', scope='spin') and mf_params[i] not in self.data_names(set='generic', scope='spin'):
1905 raise RelaxError("The parameter '%s' is unknown. It should be one of %s or %s" % (mf_params[i], self.data_names(set='params', scope='spin'), self.data_names(set='generic', scope='spin')))
1906
1907
1908 if error:
1909 mf_params[i] += '_err'
1910
1911
1912 for spin in spin_loop(spin_id):
1913 setattr(spin, mf_params[i], mf_vals[i])
1914
1915
1917 """Set all simulation selection flags.
1918
1919 @param select_sim: The selection flags.
1920 @type select_sim: bool
1921 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
1922 @type model_info: int
1923 """
1924
1925
1926 model_type = determine_model_type()
1927
1928
1929 if model_type == 'all' or model_type == 'diff':
1930 cdp.select_sim = select_sim
1931
1932
1933 else:
1934
1935 spin = return_spin_from_index(model_info)
1936
1937
1938 if not spin.select:
1939 return
1940
1941
1942 spin.select_sim = deepcopy(select_sim)
1943
1944
1946 """Function to update the other model-free parameters.
1947
1948 @param param: The name of the parameter which has been changed.
1949 @type param: str
1950 @param spin: The SpinContainer object.
1951 @type spin: SpinContainer
1952 """
1953
1954
1955 if param == 's2f':
1956
1957 if hasattr(spin, 's2s') and spin.s2s != None:
1958 spin.s2 = spin.s2f * spin.s2s
1959
1960
1961
1962 if param == 's2s':
1963
1964 if hasattr(spin, 's2f') and spin.s2f != None:
1965 spin.s2 = spin.s2f * spin.s2s
1966
1967
1969 """Initialise the Monte Carlo parameter values."""
1970
1971
1972 model_type = determine_model_type()
1973
1974
1975 param_names = self.data_names(set='params', scope='spin')
1976
1977
1978 min_names = self.data_names(set='min', scope='spin')
1979
1980
1981 if model_type == 'diff' or model_type == 'all':
1982
1983 if cdp.diff_tensor.type == 'sphere':
1984 diff_params = ['tm']
1985
1986
1987 elif cdp.diff_tensor.type == 'spheroid':
1988 diff_params = ['tm', 'Da', 'theta', 'phi']
1989
1990
1991 elif cdp.diff_tensor.type == 'ellipsoid':
1992 diff_params = ['tm', 'Da', 'Dr', 'alpha', 'beta', 'gamma']
1993
1994
1995
1996
1997
1998
1999 if model_type == 'diff' or model_type == 'all':
2000
2001 for object_name in diff_params:
2002
2003 sim_object_name = object_name + '_sim'
2004
2005
2006 if hasattr(cdp.diff_tensor, sim_object_name):
2007 raise RelaxError("Monte Carlo parameter values have already been set.")
2008
2009
2010 for object_name in min_names:
2011
2012 sim_object_name = object_name + '_sim'
2013
2014
2015 if hasattr(cdp, sim_object_name):
2016 raise RelaxError("Monte Carlo parameter values have already been set.")
2017
2018
2019 if model_type != 'diff':
2020 for spin in spin_loop():
2021
2022 if not spin.select:
2023 continue
2024
2025
2026 for object_name in param_names:
2027
2028 sim_object_name = object_name + '_sim'
2029
2030
2031 if hasattr(spin, sim_object_name):
2032 raise RelaxError("Monte Carlo parameter values have already been set.")
2033
2034
2035
2036
2037
2038
2039 for object_name in min_names:
2040
2041 if not hasattr(cdp, object_name):
2042 continue
2043
2044
2045 sim_object_name = object_name + '_sim'
2046
2047
2048 setattr(cdp, sim_object_name, [])
2049
2050
2051 sim_object = getattr(cdp, sim_object_name)
2052
2053
2054 for j in range(cdp.sim_number):
2055
2056 object = getattr(cdp, object_name)
2057
2058
2059 sim_object.append(deepcopy(object))
2060
2061
2062 if model_type == 'diff' or model_type == 'all':
2063
2064 cdp.diff_tensor.set_sim_num(cdp.sim_number)
2065
2066
2067 for object_name in diff_params:
2068 for j in range(cdp.sim_number):
2069 cdp.diff_tensor.set(param=object_name, value=deepcopy(getattr(cdp.diff_tensor, object_name)), category='sim', sim_index=j)
2070
2071
2072 if model_type != 'diff':
2073 for spin in spin_loop():
2074
2075 if not spin.select:
2076 continue
2077
2078
2079 for object_name in param_names:
2080
2081 sim_object_name = object_name + '_sim'
2082
2083
2084 setattr(spin, sim_object_name, [])
2085
2086
2087 sim_object = getattr(spin, sim_object_name)
2088
2089
2090 for j in range(cdp.sim_number):
2091
2092 sim_object.append(deepcopy(getattr(spin, object_name)))
2093
2094
2095 for object_name in min_names:
2096
2097 sim_object_name = object_name + '_sim'
2098
2099
2100 setattr(spin, sim_object_name, [])
2101
2102
2103 sim_object = getattr(spin, sim_object_name)
2104
2105
2106 for j in range(cdp.sim_number):
2107
2108 sim_object.append(deepcopy(getattr(spin, object_name)))
2109
2110
2112 """Return the simulation chi-squared values.
2113
2114 @keyword index: The optional simulation index.
2115 @type index: int
2116 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
2117 @type model_info: int
2118 @return: The list of simulation chi-squared values. If the index is supplied, only a single value will be returned.
2119 @rtype: list of float or float
2120 """
2121
2122
2123 model_type = determine_model_type()
2124
2125
2126 if model_type == 'all' or model_type == 'diff':
2127 return cdp.chi2_sim
2128
2129
2130 else:
2131
2132 spin = return_spin_from_index(model_info)
2133
2134
2135 return spin.chi2_sim
2136
2137
2139 """Return the array of simulation parameter values.
2140
2141 @param index: The index of the parameter to return the array of values for.
2142 @type index: int
2143 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
2144 @type model_info: int
2145 @return: The array of simulation parameter values.
2146 @rtype: list of float
2147 """
2148
2149
2150 inc = 0
2151
2152
2153 model_type = determine_model_type()
2154
2155
2156 param_names = self.data_names(set='params', scope='spin')
2157
2158
2159
2160
2161
2162 if model_type == 'diff' or model_type == 'all':
2163
2164 if cdp.diff_tensor.type == 'sphere':
2165
2166 if index == 0:
2167 return cdp.diff_tensor.tm_sim
2168
2169
2170 inc = inc + 1
2171
2172
2173 elif cdp.diff_tensor.type == 'spheroid':
2174
2175 if index == 0:
2176 return cdp.diff_tensor.tm_sim
2177 elif index == 1:
2178 return cdp.diff_tensor.Da_sim
2179 elif index == 2:
2180 return cdp.diff_tensor.theta_sim
2181 elif index == 3:
2182 return cdp.diff_tensor.phi_sim
2183
2184
2185 inc = inc + 4
2186
2187
2188 elif cdp.diff_tensor.type == 'ellipsoid':
2189
2190 if index == 0:
2191 return cdp.diff_tensor.tm_sim
2192 elif index == 1:
2193 return cdp.diff_tensor.Da_sim
2194 elif index == 2:
2195 return cdp.diff_tensor.Dr_sim
2196 elif index == 3:
2197 return cdp.diff_tensor.alpha_sim
2198 elif index == 4:
2199 return cdp.diff_tensor.beta_sim
2200 elif index == 5:
2201 return cdp.diff_tensor.gamma_sim
2202
2203
2204 inc = inc + 6
2205
2206
2207
2208
2209
2210 if model_type == 'all':
2211
2212 for spin in spin_loop():
2213
2214 if not spin.select:
2215 continue
2216
2217
2218 for param in param_names:
2219
2220 if index == inc:
2221 return getattr(spin, param + "_sim")
2222
2223
2224 inc = inc + 1
2225
2226
2227
2228
2229
2230 if model_type == 'mf' or model_type == 'local_tm':
2231
2232 spin = return_spin_from_index(model_info)
2233
2234
2235 if not spin.select:
2236 return
2237
2238
2239 for param in param_names:
2240
2241 if index == inc:
2242 return getattr(spin, param + "_sim")
2243
2244
2245 inc = inc + 1
2246
2247
2249 """Return the array of selected simulation flags for the spin.
2250
2251 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
2252 @type model_info: int
2253 @return: The array of selected simulation flags.
2254 @rtype: list of int
2255 """
2256
2257
2258 model_type = determine_model_type()
2259
2260
2261 if model_type == 'all' or model_type == 'diff':
2262 return cdp.select_sim
2263
2264
2265 else:
2266
2267 spin = return_spin_from_index(model_info)
2268
2269
2270 if not spin.select:
2271 return
2272
2273
2274 return spin.select_sim
2275
2276
2278 """Skip certain data.
2279
2280 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices).
2281 @type model_info: int
2282 @return: True if the data should be skipped, False otherwise.
2283 @rtype: bool
2284 """
2285
2286
2287 model_type = determine_model_type()
2288
2289
2290 if (model_type == 'mf' or model_type == 'local_tm') and not return_spin_from_index(model_info).select:
2291 return True
2292
2293
2294 return False
2295
2296
2297
2299 """Empty class to be used for data storage."""
2300